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Introduction 


The work undertaken during and this contract and its results are described in the remainder 
of the report. Fortunately many of the results from this investigation are available in the journal 
or conference proceedings literature — published, accepted for publication, or submitted for publi- 
cation. For these we simply give the reference and the abstract. The papers themselves have 
been separately delivered to NASA/GSFC. Those results that have not yet been submitted 
separately for publication are described in considerable detail. 

The accomplishments duting this contract are summarized in the following list. They 
correspond to the objectives of the revised proposal. 

[1] Analysis of the snow reflectance characteristics of the Landsat Thematic Mapper, including 
spectral suitability, dynamic range, and spectral resolution. 

[2j Development of a variety of atmospheric models for use with Landsat Thematic Mapper 
data. These include a simple but fast two-stream approximation for inhomogeneous atmo- 
spheres over irregular surfaces, and a doubling model for calculation of the angular distribu- 
tion of spectral radiance at any level in an plane-parallel atmosphere. 

(3] Incorporation of digital elevation data into the atmospheric models and into the analysis of 
the satellite data. 

[4] Textural analysis of the spatial distribution of snow cover. 
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1, Snow Reflectance from Landsat-4 Thematic Mapper 
This paper has been published. The reference citation is: 

Dozier, J., “Snow reflectance from Landsat-4 Thematic Mapper," IEEE Transac tions on Geosci- 
ence and Remote Sensing , vol. GE-22, pp. 323-328, 1984. 

Abstract. In California 75 peicent of the agricultural water supply comes from the melting 
Sierra Nevada snowpack. Basin-wide spectral albedo measurements from the Landsat-4 Thematic 
Mapper (TM) could be used to better forecase the timing of the spring runoff, because these data 
can be combined with solar radiation calculations to estimate the net radiation balance. The TM 
is better-suited for this purpose than the Multispectral Scanner because of its larger dynamic 
range. Saturation still occurs in bands 1-4, but is severe only in TMi (0,45-0,52/im). Snow 
reflectance in TM2 (0.43-0. 61pm) is typical of the visible wavelength region, where reflectance is 
almost insensitive to crystal size but sensitive to contamination. TM4 (0.78-0. 90pm) allows esti- 
mation of effective optical grain size and thereby spectral e' >nsion throughout the near-infrared. 
TM5 (1.57-1. 78pm) can discriminate clouds from snow. 




j— me— 
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2» Registering Thematic Mapper Imagery to Digital Elevation Models 
This paper has been published. The reference citation is: 

Frew, J., “Registering thematic mapper imagery to digital elevation models,” in Proceedings, 
Tenth International Symposium on Machine Processing of Remotely Sensed Data , with Spe- 
cial Emphasis on Thematic Mapper Data and Geographic Information Systems, ed. M. M. 
Klepfer and D. B. Morrison, pp. 432-435, Purdue University, West Lafayette, IN, 1984. 

Abstract. Several problems arise when attempting to register Landsat Thematic Mapper (TM) 

data to U.S. Geological Survey digital elevation models (DEMs). Chief among these are: 

• TM data are currently available only in a rotated variant of the Space Oblique Mercator 
(SOM) map projection. Geometric transforms are thus required to access TM data in the 
geodetic coordinates used by the DEMs. Due to positional errors in the TM data, these 
transforms require some sort of external control. 

• The spatial resolution of TM data exceeds that of the most commonly available DEM data. 
Oversampling DEM data to TM resolution introduces systematic noise. Common terrain 
processing algorithms (e.g. slope computation) compound this problem by acting as high- 
pass filters. 


3. Reflectance Measurements front Landsat Thematic Mapper over Rugged Terrain 
This paper has been published. The reference citation is: 

Dozier. J., “Reflectance measurements from Landsat Thematic Mapper over rugged terrain/' in 
Proceeding*, Tenth International Sympoiium on Machine Processing of Remotely Sensed 
Data, with Special Emphasis on Theme tie Mapper Data and Geographic Information Systems , 
ed. M. M. Klepfer and D. B. Morrison, pp. 230-234, Purdue University, West Lafayette, IN, 
1984. 


Abstract. Spectral albedo measurements from the Lnndsat-4/5 Thematic Mappers require that 
spacecraft upwelling radiances be corrected for atmospheric absorption and scattering and for 
local surface illumination. A two-stream model is developed, with a lower boundary condition 
that varies with incidence angle. TM data must be registered to digital terrain data. Reflectance 
from points in shadows can be used to estimate optical depth. Our primary application is deter- 
mination of the spectral albedo of snow. The TM is better-suited for this purpose than the MSS 
because of its larger dynamic range. 
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4. Automated Basfn Delineation from Digital Elevation Data 
This paper has been published. The reference citation is: 

Marks, D., J. Dozier, and J. Frew, “Automated basin delineation from digital elevation data,” 
Geo-Processing, vol. 2, pp. 290-311, 1084. 

Abstract* While digital elevation grids are now in wide use, accurate delineation of drainage 
basins from these data is difficult to efficiently automate. We present a recursive “order N ” solu- 
tion to this problem. No point in the basin is checked more than once, and no points outside the 
basin are considered. Two applications for terrain analysis and one for remote sensing are given 
to illustrate the method, using a basin with high relief in the Sierra Nevada. This technique for 
automated basin delineation will enhance the utility of digital terrain analysis for hydrologic 
modeling and remote sensing. 
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6. Field end Laboratory Measurements of Snow Liquid Water by Dilution 

This paper has been accepted for publication. The reference citation is: 

Davis, R. E., J. Dozier, E. R. LaChapelle, and R. Perla, “Field and laboratory measurements of 
snow liquid water by dilution," Water Resources Research, 1085. In press 

Abstract. Field trials of the dilution technique for measuring snow liquid water content show 
that the refined procedure is rapid and simple. Measurements of the liquid water mass fraction 
with an absolute error of ^1.5% can be obtained by one operator at a rate of 10-15 samples per 
hour, but if the water content is low, 0-2%, the relative error can be high, Electrolytic induc- 
tivity is the preferred method for measuring concentrations, using a stock solution of 0.01 N HCI. 
The recommended amount of stock solution to add is 0.5-0.8 X the mass of the snow sample. 
Extraction of the resulting mixture of stock solution and snow liquid water is best done with a 
screened pipette, instead of by decanting. 




0. Orthographic Terrain Views Using Data Derived from Digital Elevation Models 

This paper has been submitted for publication. The reference citations ere: 

Dubayah, R. O., “Orthographic terrain views using data derived from digital elevation models,'’ 
M. A- Thesis, Department of Geography, University of California, Santa Barbara, CA, 1885. 

Dubayah, R. O. and J. Dozier, “Orthographic terrain views using data derived from digital eleva- 
tion models,” Photo grammetrie Engineering and Remote Senting. (Submitted 1085) 

Abstract. A fast algorithm is present for producing three-dimensional orthographic terrain views 
using digital elevation data and co-registered imagery. These views are created using projective 
geometry and are designed for display on high resolution raster graphics devices. The algorithm’s 
effectiveness is achieved by (l) the implementation of two efficient grey-level interpolation rou- 
tines which offer the user a choice between speed and smoothness, and [2] a unique visible surface 
determination procedure based on horizon angles derived from the elevation data set. 


7. Two-Stream Method for Radiative Transfer In Inhomogeneous Atmoepheres over 
Irregular Surfaeee 

This paper has been submitted for publication. The reference citation is: 

Dosier, J. and R. F. Miliitr, "Two-stream method for radiative transfer in inhomogeneous atmo- 
spheres over irregular surf'-cs,” Journal of Geophysical Research. (Submitted 1985) 

Abstract. Two-stream a roximations for solution to the radiative transfer equation in plane- 
parallel media can be extended to inhomogeneous atmospheres over irregular surfaces. For a 
homogeneous layer the two-stream equations are solved for an irregular bounda y condition, 
which includes topographic effects and variation of reflectance with Lumination angle. Direct 
and diffuse reflectances of this layer are then used as the boundary condition for the next upward 
layer, continuing recursively to the top of the atmosphere. Accuracy of the method compares 
favorably to more precise solutions, with standard errors of 
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t. A Component Decomposition Model for Evaluating Atmospheric Effects in Remote 
Sensing 

8.1. Introduction 

Images acquired by radiometric sensors on satellites are somewhat degraded compared to 
those from lower altitude platforms because of atmospheric effects. For better derivation of sur- 
face properties and classification of ground features, it is desirable to make constructive atmos- 
pheric corrections and retrieve the ground reflectance. It is also desirable to better understand 
the relationship between the properties of the atmosphere and surface and the upwelling radiance 
at the sensor’s level. An atmospheric model with wide wavelength coverage is also useful for 
selection of optimal bands or band combinations in future remote sensing instruments. Here we 
describe an ultraviolet-through-infrared atmospheric component evaluation model for a plane- 
parallel atmosphere-earth system with arbitrarily nonuniform albedo, either Lambertian or aniso- 
tropic. Such a model can be used for testing simpler atmospheric correction models and selecting 
new wavelength bands for future sensors. 

Many researchers [Ueno et al., 1078; Otterman and Fraser, 1979; Dave, 1980; Kaufman and 
Fraser, 1984 j have noted that the radiance L received by a remote sensor is composed of three 
components: (1) directly transmitted ground signature, (2) diffusely transmitted ground radiation 
through the atmosphere, and (3) the atmospheric radiation that would occur even over a perfectly 
absorbing and non-emitting ground. Different researchers use different terminologies for these 
three components. Hereafter, we ca 1 them L f (attenuated signal), L i (diffusely transmitted 
ground radiation), and L 0 (pure atmospheric radiance). The physical meanings of these thr$c 
components are depicted in Figure 8.1. 

L = L 9 + L 4 Lq (8.1) 

Among the^e three components, L t contains useful information that we want to retrieve 
from the remotely sensed data, while the other two degrade the satellite measurement and need to 
be removed. 

Lq is usually regarded as the result o! scattering of sunlight. This is true for visible and 
shorter wavelengths, but in the more general case we can consider this as caused by both scatter- 
ing and emission. For accurate calculation of L 0t multiple scattering should be included [Ueno et 
al., 1978). 


- 10. 


The L 4 term is difficult to calculate for a iionuniform surface, because radiation reflected 
from surface areas that are not within the instantaneous field of view (IFOV) of the sensor can 
arrive at the sensor by atmospheric scattering. Ueno et al. [1078) use the mean albedo of the 
neighborhood of the target pixel to correct for individual values. Pearce (1977) takes a 
backward-tracking Monte Carlo approach to solve this problem with given ground albedo pat- 
terns. He also introduces the concept of a point spread function and points out its importance in 
retrieving ground albedo by deconvolution. Otterman and Fraser (1979) look into the significance 
of adjacency effects by a single scattering approximation. Dave (1980) uses a more sophisticated 
version of first order scattering, the “primary scattering source function 1 ’ model, to investigate the 
atmospheric eflect caused by surface inhomogeneity. Kiang (1982) assesses the importance of 
atmospheric spreading effects on Landsat Multispectral Scanner and Thematic Mapper measure- 
ments using a procedure similar to Pearce’s. Meklcr and Kaufman (1982) develop a two- 
dimensional radiative transfer model in which a one-dimension illy varying surface can be han- 
dled. Kaufman and Fraser (1984) use Pearce’s results and investigate the eflect of L 4 on 
classification accuracy. 

While all these investigations confirm the existence of neighboring area effects caused by L 4 , 
most of the approaches are forward models only, in that they calculate the upwelling radiance at 
the top of the atmosphere, given an atmospheric profile and surface albedo distribution, but they 
cannot retrieve the surface albedo distribution, given the atmospheric profile and top-of- 
atmosphcric radiance. Kiang (1982) correctly looked into the physical meaning of the atmospheric 
spread function, but he did not further investigate its possible usage in retrieving ground 
reflectance. The method of Ueno et al. [1978] is an inverse method, but without support by the 
atmospheric point spread function, it is somewhat empirical (Dave, 1980). Moreover, their atmos- 
pheric model was composed of only two layers and might be too simple. 

Pearce's [1977] model is a forward one but can be used for the inverse problem ir, some spe- 
cial cases. The wavelength range for his model is in the solar spectrum only; atmospheric emis- 
sion is not considered. The only surface type handled is Lambertian, and the point spread func- 
tion is derived only for e nadir-pointing monochromatic sensor. Diner and Martonchik [1984) 
incorporate a spatial Fourier transform method and the standard one-dimensional radiative 
transfer technique for solving three-dimensional transfer equation for a vertically inhomogene- 
ous atmosphere sitting on an inf. nogencous non-Lambert ian plane surface. So far, this method 
is a forward algorithm. 
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The development of a model for a wider wavelength range, handling both Lambertian and 
anisotropic surfaces and producing a point spread function for arbitrary sensor angles is the task 
we consider here. Given infomation about the atmospheric profile, surface reflectance can be 
retrieved from top-of-atmospheric radiance. The model therefore achieves a partial step toward 
solution of the surface remote sensing inversion problem — retrieval of surface reflection proper- 
ties from upwelling radiance measurements alone. 


8.2. Decomposition of Remotely Sensed Radiance 

To understand the physical meaning of the components of the upward radiance, let us look 
at the decomposition in a layer-structured plane-parallel atmosphert. According to the interac- 
tion principle [Grant and Hunt, 1969] radiation is additive, and the upward radiance at the top of 
a layer, bounded by upper level r, and lower level , where j — i + 1, is the sum of three com- 
ponents: reflection of top incident radiation, transmission of bottom incident radiation to the top, 
and the upward internal source emerging at the top. In vector notation, 


L f (r, ) = R(r y ,r ( - ) L*(r. ) + T{v. ,r. ) L T (r, ) + £»(r, ,r, ) 


(8.2) 


The L/s denote the vectors of downward and upward radiances in different directions at different 
levels r, and ry . R and T are reflectance and transmittance matrices, and is the upward 
internal source sector emerging at level r, . 

This equation is a statement of radiation conservation. R, T, and are uniquely deter- 
mined by the composition and status of the layer only; they are independent of the radiation imp- 
inging on the layer from outside. The L terms or the right hand side represent radiation coming 
from above or below the layer and can be caused by emission and reflection. They may include 
any interactions between the layer concerned and the adjacent layers. The decomposition of the 
transmittance term in (8.2) is a consequence of separation of direct from diffuse transmittance. 

The first and third terms on the right hand side in (8.2) result in the L 0 term. For further 
separation, we consider general expressions for L, and L d first. Then we derive equations for the 
simpler cases. For an anisotropic inhomogeneous surface, the upward direct and diffuse com- 
ponents at the sensor, which points in direction H are: 

L, = L'(r,7,n) t-'!» (8.3) 


Li =// TAo.rfl.T.r' ,H' )„• d vt dA' 

A Q 


( 8 , 1 ) 
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fl is composed of nadir angle cos^/j and azimuth 0. fl and (if are directional vectors, and </fi is 
the differential of solid angle in direction H. Therefore J dft ^ J J d d <i>. The horizontal 
position vectors f and T 9 are expre. cd by z } y and z* for the horizontal positions of 
emerging and incident radiation pencils, and dA 1 = dz f dy f is differential area. The “upward 
point-to-point bidirectional diffuse transmittance function” is 


7,(0,F,fl;r,F' .fl' )= 


d I/i' E 4 '(t,?' )) 


(8.5) 


This defines the contribution to the top-of-atmospheric radiance increment at the horizontal posi- 
tion y in direction fl made by an upwelling irradiance increment leaving from unit surface area in 
direction fl' at the horizontal position . For a plane-parallel atmosphere and surface it is 
shift invariant: it depends on the difference (r*' - F) and not on their individual values. 

The expression for L 9 (8.3) remains the same for simpler Lambertian or homogeneous sur- 
faces, so only L 4 is discussed below. For simpler cases, we first define the fallowing quantities. 
The “upward point-to-point hemispherical-directional diffuse transmittance function” from sub- 
point y 1 to sensor is 


7<(0,fl;r,F' )= J r,(0,fl;r,F' ,fl' ) a l 1 

Q 


( 8 . 6 ) 


The “total upward plane-to-poir.t hemispherical-directi anal diffuse transmittance coefficient” is 


T i (O l {l;T)= / T 4 (0,f. ;r,F' ) dA ' 

A' 


(8.7) 


For a Lambertian surface upwelling radiance at the surface is independent of viewing direc- 
tion, i.e. L f (r,F' ,fl' ) = L T (r,F' ), and (8.4) simplifies to 


Li = Li(0,fl;r,F) = / T 4 (0,fl;r,F' -F) L\t,T' ) dA 


( 8 . 8 ) 


An averaged upwelling radiance for an atmosphere of total optical thickness rat the bottom 
position F for a Lambertian surface is 


L'(r,r) 


j T d (0,fi,T,r' -F) L\r,r' ) dA 


Ti (0,fl;r) 


(8.9) 
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For a homogeneous Lambertian surface under uniform illumination ) = L '(r). 

Therefore 

L t -L'(r) ^(0,n ; r) (8.10) 

Now it is interesting to consider an opposite configuration. If a narrow beam with irradi- 
ance E 0 is incident on a unit area at the top of the atmosphere, in the same direction as the sen- 
sor viewing axis, then /i 0 ^ P and the resulting diffuse radiance at the bottom is 

Mv" & ;0,F ( {l) = # i 0 £o Tiir,? 1 ft ;0 ( F,i1) (8.11) 

and the diffuse it radiance at that location is 

E^t,?' ;0,r\fl) = / Li l {r,r' ft ,0,T ,U) ft' i W 

Q 


=p 0 E 0 fTi(T,r' ft ; 0 ,r,n) # *' dn' (8.12) 

o 

From the reciprocity principle (van de Hulst, 1980, pp. 16-18] Ti(oftj,T' ) = Ti(r ;0,f1). 
Therefore 

E^irj 1 ;0,r ,H) = p 0 #0 7^ -F) (8.13) 

The reciprocity relations among L , E , and T are presented in Figure 8.2. 

Consider the convolution expression of L d in equation (8.4). If the coordinates are chosen 
such that the zeros of T* are at the viewing axis of the sensor each time the radiance of the con- 
cerned pixel is recorded, then a new transmission function can be defined: 

T(oft,T,r ft ) = 6(r' ) *(n-H' ) + T d ( oftr,?' ft ) (8.14) 

where 6 is the delta function. Similarly, for a Lambertian surface, an integrated transmission 
function is 

T(0ftrft ) = / r(0,fl;r,r' ft ) „' dW (8.15) 

Q 

Substituting T for T d , we have 


L - L a = f T(0ft,r,r' )L'(r,r’ )dA' 

A' 


(8.16) 
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According to reciprocity again, Tiofisf ) equals the point spread function T(r,f ;0 fi' 
Knowing L and L 0 and the point spread function we can solve for L \t,7' ) by a t • 
dimensional deconvolution procedure [Andrews and Hunt, 1977] . This is the basis of applying a 
back-tracking atmospheric point spread function in restoration of remotely sensed images. 

For a homogeneous, anisotropic surface, when the pattern of L fit ) is independent of 

horizontal location , as is shown in Figure 8.3, we have the following expression for it: 

L\r,r ,fi' ) = L. r (r,?' ) „(fi' ) ( 8 . 17 ) 

where >; is an anisotropic reflectance factor. We define 


tlJOf ir.r' )= jT d (ofij,?' ,fi' ) 
Q 


d(V 


and 


(8.18) 


L » ( T ?' ) = 7 fi' )ft' d n' (8.19) 

Q 

We can substitute t dit for T d in equation (8.14), and L J ( r ,r fl ) for L\r,T' ) in equation (8.10). 
Then solve L J (t,V ) in (8.16) by deconvolution. 

For a simpler and somewhat empirical solution, go back to L d and L, . By estimating L d 
locally, the individual upward ground radiances can be solved by 

L \t,? ,H) = (L - L 0 - L d ) e r/ “ (8,20) 


L 0 and r can be accurately calculated from a one-dimensional radiative transfer model, the 
details of which are described in a following section. The value of L d can be estimated in the fol- 
lowing way, by defining 0 d =L d /L a and /? = (/, -L 0 )/L g for a homogeneous Lambertian sur- 
face: 


, _0 d (L-L 0 ) 

L i ~ a 


( 8 . 21 ) 


8.3. Aeimuthally Dependent Plane-Parallel Atmospheric Radiative Transfer 
Multiple Scattering 


Aiimuthally dependent radiance in an absorbing, emitting, and scattering layer is governed 
by the radiative transfer equation (Chandrasekhar, 1000): 

H + L C rA - Hrfl) (8.22) 

Here, the sign convention is that the downward direction is positive, r is optical depth, and 
L (r,fl) is the radiance at level r along direction fl, which is composed of zenith angle cos' 1 /! and 
azimuth The source function J is 

J (r,fl) - -£■ fP (rflfl )L (rfl) JtY + Q (r,fl) (8.23) 

The phase function I 1 ) gives the distribution pattern of single scattering at r caused 

by a pencil of radiation incident along direction fV and scattered in direction fl. The first term 
on the right hand side of (8.23) is then the total contribution made by radiation coming from all 
directions to the radiance at a particular direction The phase function is calculated by rapid 
Mic algorithms [Wiscombe, 1980], 

The Q term in (8.23) represents an internal source. By separating direct from difTuse radia- 
tion, it is convenient to consider the radiation scattered from the direct beam and the specularly 


reflected direct beam as caused by some internal “pseudo-source” [Wiscombe, 1976a], Then the 
total internal source is 

Q (r.H) = Q t (v.fl) + Q, (r,H) 4- Q, v (r,H) (8.24) 

where Q t is the thermal source and Q 9 and Q 9p are the direct and specular “pseudo-sources.” 
B[T(r)] is the Planck function. 

= (l-w)B[T(r)] (8.25) 

Q, (r,H) = P (r.titfo) < '' ,H (8.26) 

Q, P (r,H) = jL P , p (p Q ) P (r,n ; n, p ) e (8.27) 

/i 0 is the cosine of the solar zenith angle, E 0 is the solar irradiance incident on the top of the 


atmosphere (normal to the beam), p 9p is the directional specular reflectivity at the surface 
beneath the atmosphere, and r §p is the total optical thickness from top to the specular surface. 
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To resolve phase functions with strongly forward reaks using lower order polynomial 
approximation, a delta-M transformation is performed for the phase function moments, optica) 
depth, and single scattering albedo (Wiscombe, 1977). 

Interaction Principle 

One way to attack this integro-diflerential equation (8.22) uses the previously mentioned 
“interaction principle” (Grant and Hunt, 1969). In vector form, its expression for both upward 
and downward outward radiances from any arbitrary layer bounded by upper boundary r,- and 
lower boundary r y , appears as 

L*(ry ) — R(r, ,r y ) L T (r y ) + T(r y ,r, ) L^r, ) + ,r ; ) (8.28) 

L T (r< ) - R(>/ J> ) L‘(r, ) + T(r, ,r y ) L f (ry ) + S T (r, ,r y ) (8.29) 


Radiances are vectors of m X n elements on a discrete angular space composed of m zenith 
and n azimuth angles at a given optical depth: 


L"(r) = 


L{r,±n x ,<t> x ) 
L (r,±/i„<6 2 ) 


[L(r,±n m A) 


(8.30) 


1 > jii > • • • > /i m > 0 are a set of quadrature points on (0,1) and 0 < <j> { < • • • < <f> n < 2tt 
are equally spaced points in the interval 0-2 tt. The It's and T’s are reflection and transmission 
matrices, and E AT are internal source vectors. For a homogeneous thin layer, these quantities can 
be derived by some initialization scheme [Wiscombe, 1976b). 

Adding/Doubling Method 

By applying the interaction principle to two adjacent layers, the reflection and transmission 
matrices and the source vectors for the combined layer can be derived if the corresponding quanti- 
ties are known for each of these two layers [Grant and Hunt, 1969). 

Consider two adjacent layers with identical scattering properties bounded by planes at opti- 
cal depths fj, r 2 , and r 3 . By the interaction principle, we have expressions for L*(r 2 ). V(t x ), 
L^rs), and Since r u r 2 , r 3 are entirely arbitrary, we consider a single layer bounded by r x 

and r 3 , and we have new expressions for L*(r 3 ) and L T (rj). Both old and new forms for L A (r 3 ) and 
L T (r,) must be equivalent. Eliminating L 1 ^) from the first set of expressions yields the reflection 
and transmission matrices and the internal source vectors for the combined layer in terms of 
quantities for the separate layers. 
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R ( r a> r i) — Rfo.'i) + T(r,,r 8 ) [I - R(r 3 ,r 8 ) R(r,,r 8 ))-' R(r 3 ,r 8 ) T(r 8 ,r,) (8.31) 

R-(r,,r») — R(r 8 ,r 3 ) + T(r 3 ,/ 8 ) (I - R(r,,r 8 ) R(r 3 ,r 2 )]- 1 R(r,,r 8 ) T(r 8 ,r 3 ) (8.32) 

T(r s ,r,) — T(r 3 ,r 8 ) [I - R(r,,r 8 ) R(r 3 ,r 8 )]-' T(r 8l r,) (8.33) 

T(r„r s ) -* T(r,,r a ) [I - R(r 3 ,r 8 ) R(r„r 8 )]-' T(r 8 ,r 3 ) (8.31) 

“ T(r 3 ,r 8 ) (I - R(r lt r 2 ) R(r 3 ,r 8 )j 1 R(r|,r 8 ) E*(r 8 ,r 3 ) + 

T ( r s> r a) [I - R(r t ,r 8 ) R(r 3 ,r 8 )]-‘ E‘(r,,r s ) + EHrj.fj) (8.35) 

** T(r|,r 8 ) [I - R(r 3 ,r 2 ) R(ri,r 8 )| 1 R(r 3 ,r 8 ) E*(f|,r a ) + 

T(ri,r 8 ) [I - R(r 3 ,r 8 ) R(r,,r 2 )] _l E f (r 2 .r 3 ) + EVi.^a) (8.36) 


These are formulae for the “adding" method. If the two layers have identical optical thickness, 
the simpler “doubling" method results. If the initial layer is chosen such that 
fa+i ~ r i )/2 ( where N is an integer and (tv +(->■,• ) is the optical depth of the layer in the 
multi-layer system, then the reflection and transmission matrices and source vectors f r the homo- 
geneous thicker layer can be built up quickly by “doubling" N times. Note that internal sources 
are not constant with optical depth and need to be treated separately (Wiscombe, 1 076a] . 

Calculation of the Internal Radiance 

Knowing the reflection and transmission matrices and source vectors for each layer in the 
multi-layer system, we can build the internal radiance field in the atmosphere by the adding 
method. Using the formulae of the interaction principle, we have a set of simultaneous equations 
for levels 0<i <k , where k is the total number of layers in the system: 

L H r i+i) = T(r,' + i,r,- ) L l (r< ) + R(r,- ,Ty + 1 ) L T (r,' + |) + E*(r< ,r,'+i) (8.37) 

LVt ) = R( r . +i» r i ) L‘(ry ) + T(r, ,r, +1 ) L t (r < +1 ) + E T (r,' ,r f +I ) (8.38) 

The top and bottom boundary conditions that need to be satisfied are that L J (r 0 ) must be 
specified and 

L T (r* ) — Ro L‘(r* ) + t B[ T a ] + ^ « " f * /t ‘° f r (// 0 ) (8.39) 

Ro is the surface diffuse reflection matrix, ? is the emissivity vector, f r (/i 0 ) is the BRDF (bidirec- 
tional reflectance-distribution function) vector for the incident beam, and T a is the temperature 
of the surface. 
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A method to solve this set of equations is given by Grant and Hunt [1069], Its essence is 
that the radiation is additive, so that we can first consider the case in which there is no radiation 
from the bottom of the current layer to get a partial radiance, and then take in*o account the 
remaining contribution caused by upward radiance from the bottom of each layer. 


The method includes two passes. In the forward pass, the calculation starts at the top 
layer, then goes down. For each layer, only the partial radiances are calculated, which include 
the contribution made by the internal sources of the current layer plus that resulting from the 
downwelling radiation. Also, the cumulative reflection matrix and transmission-reflection matrix 
looking from the bottom are calculated for later use. This process is carried out until the ground 
surface is reached. At this point, a downward radiance over a nonreflecting, nonemitting surface 
Lb*( r * ) has been obtained. For other surfaces the interaction between atmosphere and ground is 


H(r* ) - (I - R a Ro]-‘ {W(r* ) + Ra r,(,i 0 ) + Ai( T a ) J (8.40) 


where R A is the reflection of atmosphere looking from the bottom. The downward radiance can 
then be calculated from the bottom boundary condition (8.39). 

This is followed by the backward pass carried out upwards, in which the contribution from 
the bottom of each layer is added to the previously computed partial radiances. 

The upward radiance at the top of the atmosphere is 

L'(ro) - (R(r* ,r 0 ) + T(r 0l r* ) R a (I - R(r 0 ,r* ) Rq ]" 1 T(r* ,r 0 )} L»(r 0 ) + 

S’(ro.r* ) + T (r 0) r* ) R Q (I - R(r 0 ,r* ) Rq]' 1 X 

{2W* ) + R('o.r* )\TB(T a )+ ^ e ' , ‘ / "°f r (/ 1 0 )]) (8.41) 

7T 


Fourier Transformation 

In the azimuthally dependent case all vectors are organized in lexicographic ordering 
[Andrews and Hunt, 1977]. The related square matrices are matrices with circulant blocks [Davis, 
1979]. For simplicity, we call them local circulant matrices. The results of operations of addi- 
tion, scalar multiplication and matrix multiplication of local circulant matrices are still local cir- 
culant. Moreover if the inverse of a local circulant matrix exists, it is also a local circulant 
matrix. For operations involving such matrices and vectors, the Fourier transform can be used to 
save computation time. Since the matrices are only local circulant instead of complete or block 
circulant, the Fourier transform is performed locally to take care of the azimuthal dependence, 
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and the resulting coefficients of each order of Fourier transform can also be organised as matrices. 
The sise of the resulting matrices and the length of the corresponding vectors are equal to the 
number of discretisations in the asimuth domain. 

For Fourier transforms on discrete data, the general formulae to compute the coefficients for 
sine and cosine transforms are (Scheid, 1068]: 

F„ (* +1) — A E t{j + l)sin(2fl7* /n ) (8.42) 

j-o 

Fet (* -4-1) — A E f(y + 1 )cos(2>r;A/n ) (8.43) 

/-o 

where k =*0,1,2, * • * . For even functions, such as the phase function in the present work, all 

w 

the sine coefficients are 0. The formula for the inverse Fourier transform is: 

f(j +1) == d {F e< (l) + ^i(f + l)(-l) ; + (8.44) 

f- 

2 E |F*» (s +1) coa(2irjk/n ) + F„ (k +1) sin(2jrj*/n )]} 

k — 1 

The product of the coefficients h and d in the above formulae should meet the relation 
hi =i/n, If we choose h =1, and d = 1 /ft , then for a unit matrix, the F cl ’s are all 1 for the 
diagonal subblocks and 0 elsewhere, i.e. the Fourier transform of an identity matrix is also an 
identity matrix. For computation, the Fast Fourier Transform method is used. 

With A «= 1, the Fourier transformation of a locally even circulant matrix is isomorphic. In 
other words, the forms of the original formulae remain unchanged, with the order of the matrix 
reduced. Under such transformation, the isomorphism covers the matrix operations of addition, 
scalar multiplication, multiplication, and inversion. This technique is equivalent to those used by 
Hansen and Travis [1974] and Dave and Cazdag [1970]. The computation time is reduced 
dramatically with the error introduced by the transformation of less than 10’ 7 . Note that the for- 
mula for the azimuthally averaged case is only the 0 r * order expression of the Fourier transform 
of the azimuthally dependent case. 
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8.4. Incorporation of LOWTRAN Calculations 

A a described thus far, the model is for the monochromatic case only. To make the model 
work for the atmosphere, we need to know the atmospheric optical properties. Among them the 
most important are optical thickness, r, single scattering albedo, w, and the scattering phase func- 
tions. 

The first two are related to the direct transmittance of the atmosphere. Given an atmos- 
pheric profile (temperature, pressure, water vapor density, ozone density, and the aerosol density 
a?id distribution) the LOWTRAN codes [Kneizys et al., 1983] and Mie scattering calculations give 
the atmospheric transmittance profile for wavelengths from 0.25-28.5 pm for every 20 cm" 1 
wavenumber interval. Unfortunately, the required r and w can not always be derived simply from 
the results of LOWTRAN, because a simple derivation makes the relation between the vertical 
optical thickness and slant optical length violate the cosine law that is essential for a one- 
dimensional radiative transfer model. The reason for this is that LOWTRA '' does not really give 
monochromatic transmittance but instead averaged quantities over 20 cm* 1 wavenumber intervals. 
This averaging causes violation of the Lambert-Bouguer-Betr law because of the complexity of 
molecular band absorption in longer wavelengths, even in a narrow wavenumber interval like 
20 cm" 1 . Since the total transmittance of a layer is the product of the average transmittances 
owing to molecular band absorption, molecular scattering, aerosol extinction, and molecular con- 
tinuum absorption, the problem resulting from molecular band absorption causes trouble in calcu- 
lation of the total transmittance and single scattering albedo for each layer. 

A solution to this problem is the “exponential-sum fitting” method [Wiscombe and Evans, 
1977] for radiative transmission functions calculated from LOWTRAN. For each of the three 
major absorbers (water vapor, ozone, and uniformly mixed gasses, which include C0 2 , N 2 0, CH 4 , 
CO, 0 2 , and N 2 ) the exponential-fitting is performed to get equivalent absorption coefficients k\ , 
and weights a, in this model, such that transmittance T m0 for a given absorber u is 

r m ,(«i)« (8 i5) 

i -I 

When more then one major absorber exists, the combined effect, assuming random overlap 
of absorption lines from different absorbers, is 

N . 

r m ,(u)~ n v * 

j mm \ | * 1 


(8.48) 
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Here j represent* one of N absorbers. U is a vector composed of *!,••• ,Us . For the j ik 
absorber, th, re are A/, expansion terms. Therefore, by expansion the total number of terms is 
K — A/| • • • Ms , with each having its own weight and power coefficient. The new weights and 
coefficients are 

(«") 

V - Eh, (8.48) 

for • * • ,A/y | and ro «—I, * • • % K . For each of these terms in the expansion, the mono- 

chromatic radiative transfer model can be used exactly. 

9. Model Performance 

Comparison with Oslfik and Shouman [1080] 

To verify our numerical code for a homogeneous lower boundary, we compare our results 
with those obtained by different methods. Ozi^ik and Shouman (1080) presented exact solutions, 
calculated by the Fn method, of hemispherical reflectance and transmittance values for isotropic 
incident radiation on a two-layer model with isotropic scattering properties. Stamnes and Conk- 
lin [1084] compared their discrete ordinate method with the same calculations. Now we use the 
same calculations to verify our method over a variety of single scattering albedoes and optical 
thickness combinations. We use 4, 8, and 16 discretizations in the zenith domain and analytically 
integrate over azimuth for these comparisons. In Tables 1 and 2 we show the exact solutions and 
the 4-, 8-, and 16-stream results produced by different methods for reflectance and transmittance. 
Tables 3 and 4 offer the comparison of reflectances and transmittances for the same model with 
an underlying semi-transparent specular reflecting surface. It is shown that the 8- and 16-stream 
results from our model agree with exact solutions up to 3 or 4 decimal places. For intensities 
ac UfAcy will be reduced by about one significant figure. Our results match those of Stamnes and 
Conklin (1084) to 4 decimal places in all cases but two; these minor exceptions are noted in Table 
2 . 

Comparison with LOWTRAN8 

In Table 5, we compare our results with LOWTRAN6 for the spectral interval 3-4 /im. The 
resulting upward radiances from LOWTRAN6 correspond to our results obtained when the con- 
structive contribution of atmospheric scattering is intentionally omitted. 
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Computation Spaad 

To evaluate efficiency of our code, we use a 100-layer atmospheric model with total optical 
depth of 2 for computing radiance with 10 stream*. Our model, coded in the C language, requires 
80 sec CPU time on a Digital VAX 11/780 computer operating under the UNIX operating system. 
It is difficult to compare this value with other reported times, because the run time depends on 
the computational environment. But we can at least say that this speed is comparable with those 
of Stamnes and Conklin’s [1084) discrete ordinate method, which takes 2 min 20 sec for the same 
atmosphere on the same model computer. Their model's computation time is independent of the 
total optical depth, but in remote sensing applications we are usually interested only in atmo- 
spheres with modest optical depths. Our code has the flexibility that the run time can be spent 
only on changed layers as long as we save intermediate results. 

0.1. Atmospheric Point Spread Function 

For a detailed investigation of diffusely transmitted ground radiation, the validity of a 
model using averaged ground albedo is open to challenge. Therefore, the atmospheric point 
spread function t studied. The PSF is the distribution pattern of the transmitted radiation of a 
pin-narrow beam through a degrading system. It is equivalent to the transmission in a three- 
dimensional model. 

In the present model, the atmospheric point spread function is calculated by sparse sampling 
at some specific radii and polar angles. Then the results of the sampling are smoothed by a 
least-squares fitting procedure, and a rectangular PSF is produced from the parameters describing 
the curve chosen. 

Point Spread Function Sampling Procedure 

The procedure starts by shooting photons from the receiver in a specified direction. The 
PSF is then sparsely sampled on the ground in a polar coordinate system. For each sample ele- 
ment, two distinct sampling procedures arc performed. For a Lambertian surface only the total 
number of photons that hit the clement is recorded. For an anisotropic surface, the angular dis- 
tribution of transmitted photons is recorded for each discrete direction. The angular discretiza- 
tion is exactly the same as for the radiative transfer model described in the previous section. 

Single Scattering Approximation of Point Spread Function 

For a thin atmosphere, multiple scattering is negligible, therefore a fint order scattering 
approximation is appropriate. The approximation is similar to Dave's [1980] primary scattering 
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source function model except for the following three points: (1) Dave’s approach is forward* 
tracking while this mode) is backward*tracking |Pearce, 1077]. (2) Dave's is ground albedo pat- 
tern dependent. (3) Dave does the complete integration over the entire ground surface in contrast 
to our sparse sampling approach. 

From a beam of N photons originating at the sensor, the number lost from the beam travel* 
ing between r/jt and (r-f Ar)/jj is 

AN, (r,r+Ar) - N \t -I* - e (0.40) 

Among these, the number lost by scattering is 

AN, (r,rt-Ar) AN, (r,rf Ar) w(r,r+ Ar) (0.50) 


The contribution of these photons to the ground sample element is 

AR 


A N(r, ,F,iV ) — AN, (r.r+Ar) X 

in 


/ , (t+ .ll.fi' ) exp 


r, - r- Ar/2 

..i 


(P51) 


where the exponential factor is due to the attenuation between the scattering location and the 
sampling element. AO is the solid angle increment covering the sampling element with respect to 
the point r+ Ar/2. When and Afl is almost in, special care needs to be taken to avoid 

exaggeration of the contribution because of the strong forward peak of the phase function. 

For a Lambertian surface, the contributions made by different intervals for a particular sam- 
pling element are added to get the total contribution of the beam at the location. For an aniso- 
tropic surface, the contribution from different intervals are grouped according to the angular 
discretization of the hemisphere. In this way, the point spread function for a single scattering 
approximation is calculated. 

Monte Carlo Method for Multiple Scattering 

The essence of the Monte Carlo method is that the scattering and absorption of photon bun- 
dle can be statistically simulated by a sequence of random collisions before finally the bundle is 
exhausted by absorption or escape. After collision, some portion of the photon bundle is 
absorbed, and the remaining portion rnay change the direction of motion by scattering. Each 
scattering or absorption is a random collision event, but the general trend is governed by proba- 
bility functions of the processes. 
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In this approach, the atmospheric optical parameters are the same as in the one-dimensional 
model. However, only the transmission problem is treated by Monte Carlo techniques, because 
the pure atmospheric radiance distribution has already been solved for a plane-parallel atmo- 
sphere. The problem is to find the contributions of the surface upwelling radiance to the signa- 
ture of a certain pixel; included are le directly transmitted radiance from the pixel and the 
diffuse transmission of the ground upweliing radiance. By the reciprocity principle, the pattern of 
the contributions made by the ground radiation can be mimicked by a reverse process, in which a 
beam of photons impinges at a given point at the top of the atmosphere and finally some of them 
hit the bottom and make a spread pattern on it, which is the point spread function. According to 
the interaction principle, the transmission of a layer depends only on its properties and has noth- 
ing to do with the incident radiation. The radiation interaction between layers can only change 
the amount of incident radiation, but can not change layer transmission functions. Therefore to 
calculate the transmission or point spread function, we need only consider the case in which a 
photon hits the ground once. 

In this model, the general procedures outlined by House and \vcry (1900} are followed with 
some improvements made by Pearce [1077) included, the concept of photon bundle and photon 
fraction, and the separation between a real scattering and a sampling of the contribution made by 
a scatlerer. The photon bundle concept looks at a photon as a bundle of photons and allows the 
investigator to deal with a fraction of the bundle instead of an unseparable whole photon each 
time. In this way, an absorbing atmosphere can be easily dealt with. The real scattering simu- 
lates the random walk process of a single photon in a scattering and absorbing layer. For each 
scatt ~ing event within the atmospheric layer the contribution to the point spread function is cal- 
culated. In other words, the sampling is not made when the photon hits the ground, instead it is 
made when a scattering occurs, because the diffuse radiance can be more accuiately calculated 
from the integration of the source contributions along the given path than by direct sampling. 
Such a sampling method requires many fewer incident photons. 

To mimic the random walk process, we need the distance ^raveled by the photon between 
two random collision events, the direction of each path, the portion of the photon bundle remain- 
ing after each collision, and the position of the photon in three-dimensional space. The following 
sections give the mathematical expressions of these events and quantities. Similar descriptions 
can be found in some representative papers [Cashwell and Everett, 1059; House and Avery, I960; 
Pearce, 1077). 
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Fre« path length. The distance traveled by photon bundle between collisions is called the 
free path length. Measured in the same manner as the optical depth, it is dimensionless and is 
called optical distance. The probability density function of noncollision for an optical distance / 

is 

p(/)=e-‘ (9.52) 

Then the probability that no collision occurs in the range of optical distance from 0 to / is 
i 

r, - / p(/' ) <//' -p 1 - e*' (9.53) 

0 

where r x is a uniforirly distributed random number between 0 and 1. This equation sets up a 
unique relation between a random number and an optical distance / : 

/ = -ln(l-r ,) ■* | In(l-r ,) | (9.54) 

Direction of scattering. For each scattering event, two independent angular variables 
can be obtained from the random process, the scattering angle 0 and the azimuth angle 4> that is 
measured in the plane perpendicular to the original direction 0 u <t> x . 

The scattering angle © is determined in *he following way. First define 
e 

P (6) - / p (©' )sin©' d 6' (9.55) 

o 

where 0' is a dummy variable and p(0' ) is the phase function for scattering angle 0' for the 
current scattering sublayer. Because the integration of p (©' ) over the range from 0 to jt is 2, we 
need to multiply P(0' ) by 0.5 to normalize it. Then we can relate such a normalized quantity 
to a random number r 2 to determine the scattering angle 0: 

P (0) = 2 r 2 (9.56) 

The angle d> is within the range from 0 to 2 jt; therefore 

<l> = 2nr v (9.57) 

where r 3 is another random number. Knowing the original direction and the scattering 

angle and azimuth © and 4>, the direction for the next path 0 2 ,0 2 can be determined from analytic 
geometry [Marchuk et at., 1980): 

cos 0 2 — cos0| cos© - sintf 1 sin© cos<l> 




(9.58) 
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cos^i sin^i sin© cos4> 4- cos^j sin© sin<t> 4 sintfj sin^i cos© 
cos0| cos0| sin© cos<t> - sin^j sin© sin4> 4 sin0| cos0| cos© 


(9.59) 


The portion of a photon bundle remaining. After a collision, a portion of the photon 
bundle is absorbed and the remaining portion changes direction by scattering. If the scattered 
portion does not travel horizontally, some of it may escape from the medium. Therefore, if the 
original portion before the scattering is / { and the scattering takes place at r within the medium 
of the total optical depth r*, then the remaining portion that is subject to the next scattering is 

/.“■«/. (*-« for 0 < (0.60) 

to 

/ 2 = w / , for 0 2 = 


The remaining portion / 2 will still travel within the medium. The travel distance related to ran- 
dom number rj can be determined from a transformed version of (9.54): 


Such a process is repeated until the remaining portion is too small to be of any significance. 

Distance of penetration of photon in the slab. In terms of optical depth, the penetra- 


tion is determined by 


Ar — / cos 0 2 

(9.62) 

r 2 = r, 4 Ar 

(9.6.3) 


where / is the optical distance, rj and r 2 are the optical depths for the two successive collisions, 
and Ar is the increment of optical depth between the two collisions. The height at which the 
collision occurs can be calculated from the relation between the optical depth and the height 
according to atmospheric profile. Suppose the heights for two successive collisions are h j and A 2 , 
the distance traveled between collisions is 

d = (A 2 - h ,) / cos0o (9.G4) 

Horizontal displacement. The horizontal displacement can then be calculated: 

Ax = d sin0 2 cos0 2 (9.65) 


— 11 
w. JJ 


(9.61) 


Ay = d sintf 2 sirxfo 


(9.66) 
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where Ax and Ay arc the increments in Cartesian coordinates. The horizontal location x^y 2 
can be derived given the location before traveling: 


X 2 == Ax + X | 

(9.67) 

y 3 = Ay + y, 

(9.68) 


Sampling the PSF. The sampling procedures are the same as that for the first order 
scattering except that the scattering can be of any order and can take place not only along the 
path of the the direct beam but also any place outside that particular path. For the anisotropic 
case, the direction of contribution can be over the entire upper hemisphere. 

Point Spread Function Smoothing, Curve Fitting and Interpolation 

The point spread function produced from the Monte Carlo method is not perfectly smeoth 
because of the statistical nature of the Monte Carlo procedure. Therefore, some curve fitting 
should be performed to apply the results hi radiance retrieval. 

One form that might be used is the normal distribution curve and its two-dimensional 
extension, for their wide u^e in the statistics and easy calculation of integrals over infinite range. 
However, we find that the normal curve gives low values on the outskirt of the PSF. A better 
choice is to look for the best fitting curve over a wide range of curves. A formula suitable for this 
purpose is: 


/ (*) 


i 

(A +B |* \ D ) C 


(9.69) 


where A , B , C , and D are the parameters to be determined in the fitting. This formula is 
chosen is for several reasons: (1) When the parameters A , B , C , and D are positive, the value of 
f (x) decreases with increasing absolute value of i. In other words, the curve takes a bell or 
near-bell shape. (2) When D is 2, A/B = 2(7-1. If A /B is an integer, it includes the / - 
distribution curves, often used in statistics The / -distribution curves, in turn, include normal dis- 
tribution curves when the degrees of freedom approach infinity. Therefore, if a normal distribu- 
tion curve or a / -distribution curve best fits the polar profile of PSF, then using this formula we 
can find it. (3) The integral of the volume under the bell surface for the x range O-oo (which 
represents the total point-plane directional-hemispherical transmittance function) is convergent 
when C >1. 

For an obliquely viewing sensor position, the PSF is not symmetric with respect to the verti- 
cal axis and is polar angle dependent. Therefore, it is not appropriate to use a single curve to fit 


the PSF in all the polar angles. Instead, for each polar angle, a specific curve is fitted. Then the 
volume of each piece of the PSF-bell for each polar angle, as well as the total volume of the non- 
symmetric bell curve can be found. 

Application of Point Spread Function in Image Processing 

After removing L 0t we have the following relation: 

g(z ,v) = L(i,y ) - L 0 «= £EM* -*!•*-* i) * + <« (* >V ) ( 97 °) 

This states that the ground contribution is a convolution of atmospheric PSF and the ground 
upward radiances plus a noise term e n (x t y) If the PSF h(x u y i) is known, L g (x,y) can be 
retrieved by deconvolution (Andrews and Hunt, 1977). 

The above expression is for an individual pixel. For an image, using the lexicographic form 
we can express a two-dimensional array as a vector by stacking columns for g f L and e* . Also, 
we can construct a matrix H from the PSF h such that the size of H is comparable to that of g, 
L and e n vectors. Then for the above relation over the entire image, we have 

g = HL g + e n (9.71) 

This is a system of linear equations. When the noises e n (x ,y ) are 0, the ground upward radiance 
is given by an inverse filter: 

L, = [H]' 1 g (9.72) 

Fourier transform techniques are often used in this inverse filtering (Andrews and Hunt, 1977]. In 
the current investigation, we use another technique, the “constrained least squares” algorithm 
[Hunt, 1973], to handle nonzero noise with the aid of the Fourier transform. When the inner pro- 
duct of e n , i.e. e n r e n , is estimated based on the mean and the standard deviation of the signal- 
to-noise-ratio of the sensor, this technique results in an estimation of ground upward radiance L g 
that gives the smoothest solution for given e n r e n . In other words, minimize L g r C T C L g sub- 
ject to [g - H tj T [g - Ht,j = e n T e„. The matrix C is produced from a two-dimensional 
Laplacian operator and is of the same size as matrix H. By Lagrangian method, the solution is 

L f = (H r H + tC 7- C]-‘H r g (9.73) 


7 is the Lagrangian factor that need not be solved explicitly. 
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0.2. Example - Atmospheric Effects in Lsndsst Thematic Mapper Images 

In this section we analyze the values of the three radiance components for six Thematic 
Mapper (TM) bands in visible and near-infrared regions, and we display the point spread function 
used for atmospheric correction of TM band-1 images. Finally, an image of expected ground 
upwelling radiances is retrieved from the remotely sensed TM band-1 image using that point 
s >read function. 

Three Radiance Components for a Standard Atmospheric Profile 

We chose the U.S. Standard Atmosphere [1976] with a 13-layer structure as the input atmo- 
sphere for our model because it represents an average condition for the mid-latitudes. The major 
properties and parameters of this atmosphere are shown in Table 6. The sun is assumed at the 
average sun-earth distance with the solar zenith angle at 53.7 \ The sensor is at the nadir posi- 
tion. 

Under such conditions, the three components L t , L d , and L 0 are calculated for different 
albedoes (1.0, 0.8, 0.5, 0.2, 0.0). For a given atmosphere L 0 depends on the incident solar condi- 
tion only, whereas L d and L 9 depend on the albedo. However, for a homogeneous Lambertian 
surface the ratios 0 9 =L # /L 9 and f3 d =L d jh 9 remain constant no matter what the albedo is. In 
the visible wavelength range we find that for each wavenumber interval of 300 cm -1 the relation 
that T d ( 0,fl;r) = 7rf(r,f1;0) holds quite well, because the absence of molecular absorption make 
those individual wavenumber intervals close to monochromatic cases. But in the near-infrared, 
the wide wavelength bands do not allow the monochromatic reciprocity relation to be applied, 
since the complexity of strong C0 2 and 0 2 absorption makes it completely unsuitable to approxi- 
mate such bands by monochromatic wavelengths. The difference is caused by the change of spe'e- 
tral distribution of the radiation within the wavelength interval concerned, after passing through 
the atmosphere once. For the purpose of atmospheric correction, the term 7^ (0,fi;r) instead of 
T d (r,H) is used, since that term mimics the upwelling transmission better. 

In Table 6, the terms L o, a, 0 9 , and 0 d for the 6 TM bands in the reflective solar spectrum 
are listed. These values are wavelength averaged, with the involved radiance values weighted by 
sensor response function T x and wavelength interval. 




\ Ly Tx d\ 


/T x rfX 

X 





( 





(9.74) 
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From Table 6, it is obvious that for shorter wavelengths, the terms L 0 and 0 4 have larger 
values than those in the longer wavelengths. Since 0 4 represents the total contribution of the 
diffusely transmitted ground radiation to the sensor-received radiance for a uniform Lambertian 
surface, it is an indicator of the magnitude of the adjacency effect. In TM bands 5 and 7, 0 4 and 
the ratio a*^f} 4 /0 t are email, and the adjacency effect may be neglected even for nonuniform, 
anisotropic surfaces. 

Point Spread Functions 

The shapes of point spread functions of sue bands of Thematic Mapper are shown in Figure 
0.4. The integrals of the point spread function for shorter TM bands (1-4) are close to their 
0 $ +04 sums (let 0 s ** 09 +04 and note that 0 = (L -L 0 )/L g ). This shows: (1) the results of 
Monte Carlo procedure are comparable to those from adding/doubling; (2) the reciprocity princi- 
ple holds for each narrow spectral interval. However, for some wavelengths within TM bands 5 
and 7, the integrals of point spread function even for narrow intervals differ considerably from the 
0 value when both values are comparatively small. The reason for this is that for such intervals, 
the monochromatic assumptions are no longer valid because of the complexity of molecular 
absorption. Our point spread function for each narrow interval is simulated by downward track- 
ing, starting with a smoothed solar spectrum, but the 0 values are calculated based on spectrally 
averaged upward transmissivities for each interval. The upward reflected photons have experi- 
enced longer atmospheric paths; therefore a higher portion have high penetration in the atmo- 
sphere than do the original downward spreading photons. Under such circumstances, a desirable 
point spread function cannot be obtained without renormalization. Fortunately, the effect of 
scattering at those two TM bands is negligible. A simpler atmospheric correction using L 0 and a 
will produce a good approximation for the radiance at ground level. 

Image Restoration 

Figure 9.5 shows the original images and the equivalent images restored by the “constrained 
least squares” technique. The resulting images in TM bands 1 and 2 are much sharper than the 
originals. Figure 9.6 also shows a set of restored images using a simpler technique in which the 
averaged neighborhood radiance values are used in association with the 0 4 value. A comparison 
shows that for TM band 2 the two restoration techniques give little difference. 
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Planck function at temperature T (W nr^m^sr" 1 ) 
solar constant (W m“Vtn' 1 ) 

downward diffuse irradiance at optical depth r f location t (W m^/im' 1 ) 
upward diffuse irradiance at optical depth r, location f (W nT 2 /4m~ l ) 
downward direct irradiance U optical depth r (W nrVm* 1 ) 
lexicographic noise vector 
noise term at pixel located at x f y 

BRDF (bidirectional reflectance-distribution function) for incidence 
angle cos~Vo> reflection angle cos and reflection azimuth meas- 
ured from azimuth of illumination (sr~ l ) 

BRDF vector for incidence angle cos" 1 ^ (sr 1 ) 

Fourier cosine coefficient 
Fourier sine coefficient 

lexicographic vector of radiances measured by sensor 
radiances measured by sensor 
lexicographic vector of point spread functions 
element of point spread function at location x ,y 
vector of downward radiances at level i (W m~ 2 /i nr*sr“ l ) 
vector of upward radiances at level « (VV ra'^m^sf 1 ) 

vector of downward radiances over a black and nonemitting surface 
(W nf^/im^sr"' 1 ) 

vector of expected or estimated ground radiances (VV nT 2 / 4 nT 1 sr~ l ) 
radiance at level r along direction fl (W m^jim^sr" 1 ) 
pure atmospheric radiance (VV nT 2 ^nT l sr~ l ) 

diffusely transmitted ground radiation at sensor’s level (VV nr^m^sr 1 ) 
ground upwelling radiance (W m^/im^sr" 1 ) 
attenuated signal at sensor’s level (VV m~ 2 /im~ l sr^ 1 ) 
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Averaged upward radiance from nonuniform Lambertian surface at 
ground level <W m’ 2 /im" , ar’* 1 ) 

matrix of phase functions 

phase function at optical depth r, from direction ft to fl' 

reflectance matrix 

surface diffuse reflection matrix 

reflection of atmosphere looking from the bottom 

horizontal position vector 

temperature of the surface ( *K) 

transmittance matrix 

hemispherical-directional upward transmittance from a point T at bot- 
tom to sensor 

directional-hemispherical downward transmittance from sensor to a 
ground point at location r* 

upward plane-to-point bidirectional diffuse transmittance function 

upward point-to-point hemispherical-directional diffuse transmittance 
function 

shift invariant upward point-to-point bidirectional diffuse transmittance 
function 

upward point-to-point bidirectional diffuse transmittance function 
downward point-to-point bidirectional diffuse transmittance function 

downward point-to-point directional-hemispherical diffuse transmittance 
function 

downward beam diffuse transmission at ground level 

total upward plane-to-point hemispherical-directional diffuse transmis- 
sion coefficient 

molecular transmittance for given absorber u 
molecular transmittance for a set of absorbers U 
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directional vector in three-dimensional space 
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single scattering albedo 

scattering azimuth angle measured in new coordinates 

azimuth angle, normally measured from direction of illumination 

directional specular reflectivity 

downward internal source vector (W m^/im^sr-l) 

upward internal source vector (W m^/im^sr -1 ) 

optical depth at level • , measured from top 

scattering angle between incident and scattering directions 
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sensor response function 
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Tabl« 9 4 Com pan*o n of Transmittance Result# for a TwoLayer Atmosphere above Specularly Reflecting 
Surface with Exact Calculations from Ozigik and Shouman [1080] (p %9 — specular reflectance at bot- 
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Table 9 5 Comparison of Upward Radiancs at Top of U S Standard Atmosphere at Night with Sams Cal 
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Table 9 6 Parameters Describing Atmospheric Effect on Radiances of Landsat Thematic Mapper bands, for 
U S Standard Atmosphere with 53 7 * Incident Solar Angle 
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Figure 9.2, Relation between T'b and L V (a) Definition of ,w' ). (b) Definition 

of Ti(r, ,0* ;0y,fJ). (c) By integrating T{ 0,F,f1;r,F' ,fi' ) over entire hemisphere and 
entire ground plane, T(0,i1;r) is obtained, (d) Definition of T[rfl). (a) arid (b) are a 
reciprocal pair, and (c) and (d) are another reciprocal pair. 
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Figure 0,3. An assumed anisotropic reflectance pattern in which the anisotropic reflectance factor 
depends on reflection zenith and azimuth only, and is independent of location. The gravel- 
shaped feature at each point represents a forward-peaked anisotropic refelctance pattern. 
The size of the feature indicates the magnitude of the surface albedo at that point. The 
similarity among the features shows the location independence of the anisotropic reflectance 
pattern. 
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Figure 9.5. Images before and after restoration by deconvolution using point spread function. In 
the right column are the restored images and in the left column are the original ones. The 
upper row is for the TMl images, whereas the lower row is for TM2 images The point 
spread functions are produced for the U.S. Standard Atmosphere. 
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Figure 0.5. Images restored using simpler algorithm vs. those by deconvolution procedure. Now 
the images produced by simpler algorithm using locally averaged radiance values appear in 
the right column. 







10. Texture Analysis of the Spatial Contiguity of Snow Cover 


10.1. Introduction 

Snow hydrologists have come to rely on satellite imagery for accurate estimates of snow 
cover over widespread and inaccessible areas. Prior to the advent of satellite data, snowmelt 
runoff models were based on point measurements, usually of snow water equivalent at index sites 
[Rango and Itten, 1076]. Subsequently, many studies have concluded that satellite measurements 
of sn <v-covcrcd area (SCA) have a significant statistical relationship to seasonal streamflow 
[Rango et al., 1970; Rango and Martinec, 1070; Shafer and Leaf, 1070]. This relationship has 
been used for long-term volumetric flow forecasts [Thompson, 1075; Rango et al., 1077) and for 
determining the timing of daily runoff [Martinec, 1975], 

The relationship between SCA and runoff depends on the time frame in which runoff is 
viewed. For daily predictions during the ablation period, SCA is directly related to runoff. The 
Martinec runoff model that has proven useful for many mountain basins [Rango and Martinec, 
1981] expresses the relationship 

Qt = t(E X SCA X (1-*) + i) (10.1) 

Here Q t is runoff at time t , e is a runofT coefficient, E is energy, usually in degree-days, SCA is 
snow-covered area, and k is the recession coefficient. 

In the longer term, it is known that the rate of snowline retreat is inversely related to snow 
water equivalent and to runoff [Rango and Itten, 1976]. Furthermore, the relationship between 
SCA and depth is quite variable: snow covering the same areal extent can vary 200 percent in 
depth (Martinec, 1980]. In fact, Martinec found that SCA is better related to the ratio of the 
current water equivalent to maximum seasonal water equivalent. This agrees with Thompson’s 
[1975] earlier finding that SCA is more strongly related to the percent of total seasonal runoff than 
to runoff itself. Thus the behavior of the SCA parameter in runoff models beyond those for 
short-term forecasts becomes complex; it is best represented in a form differentiated with respect 
to total accumulation [Martinec, 1980]. Alternatively, for long-term forecasting, the more direct 
relationship between snow depth and total water volume may be preferable. From the analysis of 
a large random field sample, Adams and Roulet [1982] found a broad similarity in patterns of 
depth and water equivalent both in terms of quantity and distribution. They suggest that depth 
may be a good indicator of water equivalent and therefor j of runoff. 
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Snow depth cannot be directly determined from a single satellite image but the pattern of 
snow distribution at various overall depths is readily apparent. In fact, expert snow observers 
have used the appearance of surface features to estimate depth during aerial overflights [Rango 
and Itten, 1970] . The notion that snow depth can be inferred from snow cover patterns implies 
that accumulation patterns are stable over time due to the control of underlying physical vari- 
ables of terrain and exposure. During the accumulation period snow cover patterns are event 
related and result from irregular deposition influenced by elevation, wind and local topography. 
On the other hand, overall patterns that form during the snowmelt period are quite predictable 
since melt rates are strongly controlled by altitude and exposure. Palmer [1981] found over a 
three year period in the Rio Grande watershed of Colorado that snowline recession patterns were 
repeated; Lichtenegger and Seidel [1981] reviewed images of the Dischma valley in the Swiss Alps 
over an eight year period and concluded that a typical snow cover pattern forms each year during 
melt season. Moravec and Danielson [1979] and Martinec [1980] have also reported that yearly 
repeating contour patterns of snow-covered regions occur during the ablation period. 

The analysis of snow cover patterns has generally been conducted as part of research into 
mesoscale (lOO-lOOOm) areal differentiation of snow cover such as those based on identifiable 
landscape units [Adams and Roulet, 1982] or hydrologic response units (Thomsen, 1980], The 
purpose of these investigations is to develop regional generalizations about snow conditions from 
sites stratified by similar combinations of environmental variables. By inversion, it can be 
reasoned that the snow patterns themselves are meaningful expressions of the sum effect of the 
controlling variables. As Palmer [1981] points out, the position of the snowline acts as a natural 
integrator of the long-term effects of snow accumulation, slope, aspect, temperature, radiation, 
and wind. In a one-dimensional approach to the problem of quantifying snow patterns, Palmer 
developed regression relationships between percent snow cover and snowline elevation along a 
series of index baselines for the purpose of predicting SCA for an entire basin especially at times 
when it is partially obscured by clouds. This method requires that the network of baselines in a 
basin include all areas of significant snow cover but ignore detached patches of snow. It would 
stem more appropriate to use a two-dimensional characterization of the spatial contiguity of snow 
cover to predict snow-covered area. 

Objectives 

In this study it is proposed that two-dimensional descriptions of snow cover obtained by 
means of texture analysis, a set of statistical pattern recognition techniques, serve as predictor 
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variables in a linear relationship for predicting snow-covered area. This is undertaken as a prel- 
iminary step to assess the feasibility of predicting in the same manner another hydrologic parame- 
ter, snow depth, which may have a more direct and consistent relationship to the total volume of 
water stored in a melting snowpack. Operationally, the estimation of snow coverage itself on the 
basis of two-dimensional statistics may be useful along transition zones of the pack and during 
melt season when snow cover is highly dissected and difficult to inventory by manual or digital 
means. At this time, when short-term forecasting occurs the SCA variable is a direct, useful pred- 
ictor of daily runoff. Snow depth data at the scale an 1 extent necessary to conduct pattern 
analysis were not available for this study, however it may be possible in future studies to photo- 
grammetrically determine snow depth at a scale appropriate to the analysis [Cooper, 1965; Rawls 
and Jackson, 1679]. 

As an initial stage in the investigation, the effect of sensor resolution on detectable “tex- 
ture” is studied to determine whether the large improvement in spatial resolution provided by the 
Landsat Thematic Mapper (TM) sensor over existing Multispectral Scanners (MSS) translates into 
equally improved spatial information when analyzed using standard texture analysis methods. 

Study Design 

Digital images from the LandsaM satellite were available at 30 meter resolution (TM) and 
at 80 meter (MSS) resolution. Four matched sets of TM and MSS subimages (Figures 10.1 and 
10.2) were selected and registered. For one set (image A), texture statistics were calculated over 
the entire image in order to closely investigate the behavior of the statistics at both resolutions. 
In the next step, two sets of texture features were calculated from windowed samples over each of 
the four image pairs. The relative distance between sample texture features was assessed by three 
different metrics. In the last stage, binary classifications of snow were made at both resolutions 
for image pair A. From these SCA was calculated by window and regressed against the two sets 
of sample texture features calculated above. Model efficiencies were calculated both internally 
using a jackknife regression technique and through time by cross-prediction between images of the 
same site having undergone significant snow recession. 

10.2. Texture Analysis 

Satellite images are two-dimensional projections of the three-dimensional landscape below. 
Frequently, such images ar* used to supply point data about scalar quantities like brightness, 
temperature and elevation or used in combination with other images to provide vector 
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information such a a color [Ahuja and Sch&chter, 1981). In doing so, the higher order relationships 
inherent in the image are effectively ignored. However, in one area of image processing, texture 
analysis, the study of spatial relationships is well developed. In general, the aim of texture 
analysis is to apply pattern recognition techniques to an image in order to segment and discrim- 
inate between scene regions or to aid in classification of cover types; the task is to quantify an 
invariant, non-labile characteristic of a scene object. Characterization of changing spatial pat- 
terns has not been well explored. The objectives of the present study are to use texture analysis 
in the traditional sense to quantify resolution-dependent differences in texture and to explore a 
new possibility of using extracted texture features as meaningful parameters in a functional rela- 
tionship for the prediction of a physical variable. 

Texture Analysis Methods 

Texture analysis, the image processing term for pattern analysis, originates from empirical 
efforts to recognize and duplicate the elusive perceptual concept of texture. Visual analogies have 
held sway so long in this field that only recently have formal image models emerged on the level 
of abstraction found in other spatial disciplines [e.g. Pielou, 1977). Within the sizable battery of 
empirically developed methods [reviewed by Haralick, 1979], no single approach has proven to be 
universal, in large part because the visual hierarchies involved in perceiving spatial structure work 
in a complex manner not easily duplicated by simple methods [Julesz, 1975). As Haralick [1979] 
has noted, the organization of tonal primitives or local regions can be viewed as structural, proba- 
bilistic, or functional depending on relative resolution. Whether stochastic pixel-based models or 
deterministic region-based models are the most suitable texture descriptors depends on the coarse- 
ness, homogeneity and periodicity of the texture. 

Statistical approaches range from simple first-order measures like grey tone differences and 
run lengths [Galloway, 1975) to more complex joint and conditional second-order co-occurrences 
[Haralick et al., 1973], One-dimensional autoregressive models [McCormick and Jayaramamurthy, 
1974; de Souza, 1982] are only partially successful at describing spatial patterns while two- 
dimensional autoregression [Tou, 1980) becomes a complex task. 

Image patterns can be analyzed in terms of spatial frequency but Fourier analysis has had 
limited application to texture analysis. The Fourier transform must be computed over large win- 
dows and comparison of power spectra between different sized regions is difficult [Chen, 1979). 
More importantly, local information is scattered in the frequency domain so that similar peaks 
may be caused by a nearly periodic texture or a single strong edge [Nevatia, 1983). Recently 
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Clicn 11082) and Jernig&n and D’Astous (1081) have successfully developed local and size invariant 
texture features based on the Fourier transform that overcome some limitations of the method, 
findings from perceptual experiments [Pratt et al, 1078; Julesz and Caelli, 1979] have cast doubt 
on the efficacy of Fourier analysis for texture discrimination. Patterns having identical power 
spectra and thus identical autocorrelation functions, can be discriminated effortlessly by eye. 

Translated into the spatial domain, Fourier analysis is simulated by a series of convolutions 
[Faugeras, 1078; Laws, 1080]. Convolution masks which enhance high frequency information act 
as edge detectors that approximate mathematical gradient operators [Ballard and Brown, 1982]. 
Once obtained, the edge structure of an image can be reported simply in spatial averages or used 
to form high level primal sketches [Marr, 1982). Related to edge analysis arc methods that quan- 
tify local maxima and minima by row [Mitchell et al., 1977] or which construct more complex 
relational trees of one-dimensional intensity profiles expressed a.«> nested or concatenated peaks 
[Ehrich and Foith, 1978). 

When the elements of a texture become much larger than the resolution cell of an image, 
pixel-based stochastic models break down and are supplanted by structural methods which iden- 
tify primitives, measure their attributes and determine their spatial relationships [Wang et al., 
1981; Matsuyama et al., 1982; Tomita et al., 1982]. In highly regular patterns, primitives can be 
described syntactically using tree grammars [Lu and Fu, 1979]. 

There have been few studies undertaken to rigorously compare texture analysis methods. 
Frequently cited works by Weszka et al. [1976] and Conners and Harlow [1980] have led to the 
widespread use of Haralick’s second-order statistical features: the moments of the grey-level co- 
occurrence matrix (GLCM). Indeed, co-occurrence statistics have been very useful for image seg- 
mentation [Chen and Pavlidis, 1979; Conners et al., 1984] and for image classification [Hallada et 
al., 1982; Vickers and Modestino, 1982; Holmes et al., 1981], Julesz’s [1975] finding that human 
texture discrimination operates finding that human texture discrimination operates at the level of 
second order relationships has lent such support for the GLCM approach that less costly first 
order methods reported to perform equally well for classification purposes [Weszka et al., 1976; 
Mitchell and Carlton, 1978; Pietikainen et al., 1983] are not implemented as often. Because the 
GLCM serves as the standard of comparison for testing the performance of texture analysis 
methods, it was chosen for use in this study along with a newly reported local method: Laws’ 
[1980] texture energy measures. 
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Gr«y-Uv«l Co-occurrancc Matrix 

The grey level co-occurrence matrix (GLCM) ia an array of joint frequencies whose dimen- 
sion is equal to the number of grey levels in the image. Each entry in the matrix is the frequency 
with which brightness i co-occurs with brightness j when separated by distance d in the direction 
#. Frequencies, often normalized to probabilities, are reported in both directions for a joint pixel 
pair making the matrix symmetric. As a first step towards data reduction, the number of grey 
levels is decreased to 04 or less, and the four directional matrices can be averaged into one. To 
compress the data further, several statistics are calculated that express either the distribution of 
matrix values around the main diagonal or the degree of correlation between matrix rows and 
columns. Seven statistics proposed by Haralick et al. [1073] are given in the Appendix; these 
include energy, correlation, homogeneity, entropy, inertia and information correlations 1 and 2. 
Energy and homogeneity are measures that emphasize low contrast transitions; entropy and iner- 
tia increase with texture coarseness. Correlation statistics measure the degree of association 
between marginal and total values expressed either as frequencies or entropies. 


An elegant solution to the problem of choosing an optimal combination of dist ince and 
orientation to best describe the structure in a texture was proposed by Zucker and ferzopoulos 
[1980] . They developed a chi-square statistic based on maximum likelihood estimates of the mar- 
ginal matrix probabilities to test the independence of rows and columns. The unnormalized co- 
occurrence matrix is thus viewed as a contingency table in which intensity pairs are samples 
obtained from a two-dimensional random process. Notationally: 
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N is the total number of samples. Degrees of freedom u — (hi -l)(n -1), i,y is the co-occurrence 
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In the present study of snow texture, an automated system was developed to calculate the follow- 
ing for each subsample of an image* iour unnormalized directional matrices, the chi-square value 
for each, the normalized matrix for the maximum chi-square angle, and the seven co-occurrence 
statistics. Only the final statistics from the most structured matrix were concatenated into an 
output file. This process was repeated for each of four distances (1,3,5,10 pixels) and four quanti- 
zations (8,16,32,64 grey-levels) for image pair A. 


- 52- 


Tixturt Energy Measures 

Despite the shortcomings of Fourier analysis outlined above, texture analysis methods that 
characterize sputial frequency are very useful if implemented locally in the spatial domain. This 
is accomplished by a sequence of boxflltcr applications which are faster and simpler than a single 
convolution using the fast Fourier transform [McDonnell, 1081]. In his dissertation, Laws [1080] 
derived a series of one* dimensional operations of certer-weighted local averaging, symmetric first 
differencing (edge detect j«m) and second differencing (spot detection) [Pietikainen et al., 1083], 
When convolved together these vectors form nine 3x3 masks (see Appendix) some of which are 
recognizable as standard gradient operators like the pair of vertical and horizontal Sobel operators 
(/•7 t , fU % ), and the Laplacian second difference operator (/t7«). Note that all but /i/ 0 , the low 
pass smoothing filter, are zero-sum filters. 

Texture features are obtained from each of the nine separately convolved images by calcu- 
lating local statistics such as the sum, the mean or the standard deviation over small windows. 
McDonnell [1081] along with Laws have found that the variance or standard deviation of filtered 
windows are very powerful measures of image texture. In a zero mean field produced by convolu- 
tion with a zero-sum mask, variance is the average of the squared values which makes it a meas- 
ure of total energy within a window. Laws claims that the average absolute valu»; is a fast 
approximation to the standard deviation; he refers to both the average and the standard deviation 
as measures of texture energy. Pietikainen et al. [1983] tested two other texture energy features, 
the sum of the absolute values and the maximum value within a window, and found that the local 
maxima performed just as well as the sum. In this present study three features were compared: 
the sum, average and the standard deviation of values within windows sized 16 x 16 on the MSS 
image and 32 x 32 on the TM image. Following Laws* convention, these features are referred to 
as SUM, AVG and SD. As was done for GLCM analysis, an automated system was developed to 
cycle through all nine filters, convolve the image, compute the local statistics by window and con- 
catenate them into an output matrix. Window size was based on Laws’ finding that classification 
accuracies were nearly perfect using 32 x 32 window but dropped rapidly below 15 x 15. Accord- 
ing to Hallada et al. [1982], this sample size is also adequate for co-occurrence analysis. They 
found that class separability increased logarithmically and then leveled off as window size 
increased from 3 x 3 to 13 x 13. 

Texture energy measures are distinguished from Fourier methods by their local nature. 
Phase relationships within each window are measured without reference to a global origin [Laws, 
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1080], whereas Fourier frequency components contain global information from across an entire 
image at the neglect of local information [Jernigan and D’Astous, 1084). In addition, gradient 
filters can be tailored to various textures by increasing the differencing distance to overcome noise 
but keeping it small enough so that local gradient remains a good representation of local changes 
[Ballard and Brown, 1082]. 

Parallel work by Faugeras and Pratt [1080] lends support to the Laws energy approach. 
Because the autocorrelation function has proven insufficient for texture discrimination (see above) 
these authors sought ways of characterizing the decorrelated texture field which would yield useful 
texture measures. Decorrelation can be accomplished by a whitening operator based on adjacent 
row and column correlations; if correlations are perfect this operator becomes the Laplacian 
operator. Alternatively, gradient operators like the Sobel filter can replace tbr whitening opera- 
tor. Note that these are three of the nine Laws convolution masks. 

Using a distance metric criterion, Faugeras and Pratt [1080] found that the first four 
moments of the first-order histogram of the decorrelated field provided good separability between 
similar natural textures. The first two histogram moments of the decorrelated images are exactly 
equivalent to the average and standard deviation of texture energy planes convolved with the 
same gradient operator. In accord with perceptual findings, inclusion of ohape measurements 
taken from the autocorrelation function improved separability but alone were weak discrimina- 
tors. The Sobel operator, a directional filter that does not zero out the mean or create unit vari- 
ance, gave the best separability while the non-direntional Laplacian was the worst. This implies 
that Law’s choice of average and standard deviation features is well founded since all but the 
smoothing and Laplacian filters are non-symmetrical. 

Pietikainen et al. [1983] have confirmed that local statistics of convolved images yield better 
classification results than co-occurrence statistics. In the following analysis of snow cover patterns 
the two methods are compared for relative powers of separability and utility in characterizing 
spatial distribution for the purpose of predicting area. 

10.3. Data Processing 

Registration and Sampling 

The imagery used in this study was taken by Landsat-4 Thematic Mapper (TM) and Mul- 
tispectral Scanner (MSS) sensors. TM imagery was available for two dates, December 10, 1982 
and January 18, 1983, in which images overlapped along adjacent paths; identical MSS imagery 


was available only for the December date. The images are centered on the Kern and Kings river 
basins in the southern Sierra Nevada which are adjacent watersheds on the order of 4,500 square 
kilometers (see Rango et al. (1979) for a complete geographic description). Four subimages of 
varying textural complexity sized 256 pixels in dimension were selected from the December TM 
image. Corresponding MSS subimages, 128 pixels in size, were located and registered to the TM 
sites. Registration was a simple matter of enlarging the MSS subimage two-fold and translating 
the image to line up with the TM subimage. Resampling was unnecessary because geometric 
rectification pe nrmed by the NASA Goddard LAS system left MSS resolution almost exactly half 
the TM resolution. This level of registration accuracy was adequate for comparison of textural 
differences between resolutions. For the second stage of regression analysis a single site in the 
Kern basin was selected from both December and January TM images that showed evidence of 
substantial snow recession between scenes. These subimages were also co-registered using simple 
translation without initial resampling. 

The texture study sites are about 65km 2 , comparable in size to several small experimental 
watersheds [e.g. Rango and Martinec, 1979). For an initial comparison of the behavior of texture 
features at the two resolutions, co-occurrence matrices weie calculated over the entire scene. For 
the purpose of separability measurements and regression analysis image pairs were subsampled 
using 32X32 sized windows for the TM image and 16X16 sized windows for the MSS image. 
This non-overlapping sampling strategy yielded 64 samples per subimage each covering about one 
square kilometer. At this scale, the analysis remains within the realm of mesoscale studies and is 
equivalent in scale to NOAA AVHRR imagery at nadir. 

Background Effects 

Much of the experimental work done in texture analysis has been carried out on homogene- 
ous texture fields which are assumed to be consistently specified by either parametric or deter- 
ministic models derived solely from the relationships of the texture primitives. Texture analysis 
of natural terrain must take into account external variables such as topography and vegetation 
that act as forcing functions on the pattern of surface cover. Shadows, topography and plant 
cover become part of a scene spe** dc textural characterization of the overlying snow cover. 

Topographic effect and shadowing can be reduced if digital elevation d*ta are available, 
making it possible to map radiance values into a a synthetic brightness image using lambertian or 
non-lambertian models of surface reflectance [Justice et al., 1981], Digital elevation data are 
available for the southern Sierra Nevada only at 90 meter resolution and according to Seidel et al. 
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[1983] elevation data on a scale greater than the Landsat sampling distance are inadequate for 
producing a satisfactory synthetic brightness image. Band ratioing, a method for canceling out 
multiplicative effects of topography was rejected because of its tendency to enhance differential 
noise between bands giving rise to spurious or confounding signals of high frequency texture. 
Without the possibility of digital terrain correction, images were selected in which shadowed areas 
were a small proportion of the image. This unfortunately limited the method of analysis to 
larger, open areas of mountain basins which may not be truly representative sites. Since all com- 
parisons in this study were scene specific, it was assumed that the texture signal caused by under- 
lying factors would hold constant between resolutions and between dates. More importantly, the 
chosen texture analysis methods, co-occurrence and energy statistics, are sensitive measures of 
contrast and of edge structure and should thus reflect the distribution of very bright, high- 
contrast snow patches rather than dark, low contrast background features. To emphasize con- 
trast and edge detad all texture features were derived from visible bands TM 2 and MSS 4 in 
which snow is very bright. 

Pre-Processing 

Texture measures, like co-occurrence statistics that are based on grey-level transitions, are 
sensitive to shifts in overall scene brightness or contrast. To standardize images so that mono- 
tonic changes in illumination are not used to discriminate textures, the first order grey-level distri- 
butions of all the textures were normalized to uniform distributions. At the same time, images 
used for GLCM analysis were reduced in quantization to make the co-occurrence matrices reason- 
ably sized. Histogram equalization was carried out using a procedure outlined by Pratt [1978]. 
This process can be considered a monotonic point transformation in which the input cumulative 
probabilities are equal to the output cumulative probabilities for a given input index. The histo- 
gram equalization function is expressed: 

9 = [ 9 “ 9 min] if ) + 9 min (10.3) 

Here g is the output grey value, 0 min is the minimum output grey value, is the maximum 
output grey value, and Pj (/ ) the cumulative distribution function of the input variable / . 

Note that the output number of grey levels is controlled by the g - g m \ n range, so that 
images are simultaneously equalized and reduced in quantization by a single transformation. 


10.4. Snow Claulfication 


Spectral Characteristics 

Snow has a very distinct spectral signature. It Is extremely bright in the visible bands, fre- 
quently saturating sensors calibrated for vegetation reflectances and at the same time is quite 
dark in the shortwave infrared bands like TM bands 5 and 7 [Dozier, 1984] Few scene elements 
are confused with snow cover except for white clouds that often are indistinguishable in the visi- 
ble and near infrared spectral range. With the advent of TM shortwave IR data (1.57-1. 78/im 
and 2.10-2.35/im) discrimination between the two classes has become possible because clouds are 
significantly brighter than underlying snow in these bands [Dozier, 1984]. This suggests that 
given visible and shortwave IR data, satisfactory snow classification could be achieved with only 
two spectral bands. 

A successful two band snow classification using MSS visible and near infrared data is 
already in use; Haefner [1979] found that snow cover could be classified into three found that 
snow cover could be classified into three categories, snow-free, transitional and snow-covered, 
using only MSS bands 5 and 7. Although it was necessary to further subdivide classes during 
training site selection, the visible and near infrared bands were sufficient for discriminating the 
three snow categories except in the presence of concrete, white rocks or snow under dense coni- 
ferous forest. (Mixed classes of snow under forest canopy can be discriminated with the addition 
of MSS band 4). For complex classifications which distinguish various stages of snow metamor- 
phism all four MSS bands have been used to identify up to ten classes of snow and seven classes 
of ice [Thomas et al., 1979]. 

Classification Approach 

Considering the superior resolution and spectral discrimination of the TM sensor, it was 
presumed that a two band approach using TM bands 2 and 5 would be an improvement over MSS 
two band methods and sufficient for a binary classification of snow-covered and snow-free areas. 
In the visible range TM band 2 was selected because of its larger dynamic range and therefore 
lower tendency to saturate over snow compared to TM band 1 and yet remain relatively insensi- 
tive to metamorphic changes in grain size compared to TM bands 3 and 4 [Dozier, 1984]. TM 
band 5 had several reasons to recommend its use; besides cloud discrimination properties, this 
band is generally a high information channel. Price [1984] found that the TM band 5 information 
rate expressed in bits/pixel is higher than shorter wavelengths and is also relatively uncorrelated 
with visible and near infrared bands. 
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Because snow is spectrally distinct, very bright and relatively homogeneous it is amenable to 
thresholding approaches to classification. In fact, Rango and Itten [1076] found that snow 
classification results differed little between histogram parallelepiped and maximum likelihood 
classifiers. For this study, rather than use a parallelepiped scheme, methods were investigated to 
find a single data plane for thresholding which combined critical spectral properties of visible and 
shortwave infrared wavelengths and also reduced variance. 

Table 10.1 gives summary statistics for a single snow cover training site from TM image A 
(n=*=315) using two approaches, ratioing and principal component analysis. Ratioing the two 
adjacent bands, TM 2 and TM 3, reduced variability and range when compared to TM band 2 
alone. Presumably, variation due to topography was reduced by canceling out multiplicative 
effects. This reduction is greatest for adjacent bands in which surface reflectance ranges are simi- 
lar [Holben and Justice, 1981]. Although the narrower threshold of the ratioed class was an 
improvement, it did not include important shortwave IR spectral information. A ratio of TM 
bands 2 and 5 reduced the threshold range but not the variability with respect to the single visi- 
ble band. 

On the other hand, snow patches were clearly discernible on principal component images. 
As Figure 10.3 shows there is virtually no difference between the first principal component using 
TM bands 2, 3, 4, 5 and 7 and the first component using bands 2 and 5 alone (see discussion below). 
The two band component was selected for thresholding since it was computationally less costly 
and because it had a lower coefficient of variation than either TM band 2, the TM2/TM5 ratio or 
the five band principal component image. A single threshold range applied to the two band prin- 
cipal component image yielded a satisfactory classified image (Figures 10.4 and 10.5). 

Classification Results 

The spectral and spatial advantages afforded by the TM sensor are evident from comparing 
TM and MSS snow classifications for the .’ame scene (Figure 10.4). The MSS scene was also 
classified by thresholding a single data plane, the first principal component of bands 4,5,6 and 7, 
which was selected by the same process of comparing training site statistics. Percent snow cover 
in the MSS subimage was 30% higher than in the matched TM subimage. This discrepancy 
partly stems from the poorer spectral discrimination provided by the MSS data and the inexact, 
subjective nature of thresholding, but its major cause is the far coarser resolution that blurs tran- 
sition zones and leads to systematic overclassification. At the same time, isolated, small groups of 
snow-covered pixels identified on the TM image were omitted from the MSS classification. 
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Principal Component Analysis 

Principal component analysis is a linear orthogonalizing transformation that yields a vector 
of coefficients, the eigenvector, used in linear combination with the variables of the input vector 

f; 

to align the variables along an axis of maximum variation that is statistically uncorrelated and 
geometrically orthogonal to Sc seeding components [Cooley and Lohncs, 1971). A principal com- 
ponent image is obtained from the original image, g , having p spectral bands by the transforma- 
tion [Moik, 1980]: 

g ' T (g - m) (10.4) 

Here g e is the principal component image, T = pip is the matrix whose rows are the normalized 
eigenvectors of the spectral covariance matrix C of g , and m is the mean vector of the p spectral 
bands. 

The eigenvalues \ p and the eigenvectors t p of C are obtained by solving the equation 

C f p = Xp t p (10.5) 

The principal component transformation isolates non-random information from noise while also 
decorrelating the transformation axes to eliminate redundancy [Anuta et al., 1984], The resulting 
scalar eigenvalues and set of eigenvectors can be interpreted directly for some insight into the 
sources of variation. The ratio of each eigenvalue to the sum of all eigenvalues gives the percent 
of total variance explained by the corresponding eigenvector. The eigenvectors are comprised of 
coefficients or loadings that correspond to the cosine of the angular distance through which each 
input band must be rotated to be aligned with the principal axis of variation. The larger 
coefficients represent smaller angular distances and thus greater influence on the component 
[Anuta et al., 1984], The loadings can be viewed as weights corresponding to the relative contri- 
[ bution made by each band. 

Component loadings for TM and MSS subimages are given in Table 10.2. Clearly, the visi- 
ble bands dominate the first principal component derived from five TM bands; visible and near IR 
; bands all carry about equal weight while the shortwave IR bands contribute far less. Only in the 

second component does the near IR band behave independently of the visible bands. The load- 
ings of the first principal component using only two bands reflects this same pattern and provides 
the same level of explained variance as the five band first component. Apparently redundancy 
among the five TM bands with respect to the first axis of maximum variation was effectively elim- 
inated by using only two bands. 
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Loadings for the first component derived from four MSS bands arc not clearly divided into 
two classes. Visible bands show a pattern of increasing relative contribution into the IR which 
then tapers off again; MSS bands 4 and 7 alone were unable to represent this combination of spec- 
tral information. In reviewing the component images, it is apparent that TM snow classification 
was achieved using essentially a reduced-variance visible image augmented slightly by shortwave 
IR information whereas MSS classification relied on a nearly equal mix of visible and near IR 
bands. 

10.6* Methods 

Distance Measures 

In the language of statistical pattern recognition, the texture statistics used in this study are 
features which detect sufficient statistical (non-random) variability between patterns to allow 
classification. The task in feature selection is to find an evaluation function that will assess how 
well a set of features discriminates between classes. Generally, there are three types of evaluation 
rules [Ben-Bassat, 1980]: information measures (uncertainty), distance measures (separability), and 
dependence measures (association). Each of these measures distribute objects into feature space 
which can be divided into classes by a discriminant function. Alternatively, each measure can 
stand as a figure of merit such that a large measurement difference implies low classification error. 
Evaluation by figure of merit rather than by classification has the advantage of being independent 
of any particular discriminant function and may additionally include error analysis [Faugeras and 
Pratt, 1980]. Distance metrics were deemed most appropriate to the aim of quantifying 
resolution-dependent textural differences rather than for discriminating between them. 

Distance measures used in pattern recognition for statistical evaluation of separability 
operate on sample pools. Accordingly, the large TM and MSS images were divided into 64 square 
samples covering comparable areas. Distances were calculated between these samples and then 
reported in sum or average. 

Three distance measures were chosen to evaluate the separability of TM and MSS texture 
features. Initially Euclidean distance was calculated between sets of statistics that were first re- 
scaled between zero and one in order to preserve a consistent metric for variables originally meas- 
ured on different scales. The Euclidean distance between columns of the TM and MSS feature 
matrices was calculated as follows: 

<*(;.*)= E '/[*(*,) )j 4 
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The Euclidean distance between image pairs was the sum total of distances between all texture 
features (matrix columns). Because it was necessary to re-scale the data, this metric accommo- 
dated comparisons between unnormalized variables such as sums and counts. 

To be able to easily compare the degree of similarity among data sets, the average Gower 
Similarity Coefficient [Gower, 1071] was calculated between each feature k: 
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R (it) is the range of a given texture feature k over both sets of data. The similarity between 
image pairs was the overall average similarity between texture features. Like Euclidean distance, 
the similarity coefficient is sensitive to magnitude. When applied to inherently normalized 
features such as the standard deviation or the average, the range is sufficiently comparable 
between data S 2 ta for the coefficient to work well. Unsealed data sets ranging widely in value will 
have low similarity though correlation between the two may be high. 


The last distance metric considered, Bhattacharyya distance ,is a more sophisticated meas- 
ure theoretically based on a scalar function of the probability densities of the two feature sets 
[Faugeras and Pratt, 1980]: For Gaussian densities Bhattacharyya distance is calculated [Davis, 
1981]: 
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Here ti t - is the mean vector for class « , E, is the covariance matrix for class t , and | E, | is the 
determinant of E,- . 

By taking variability into account, this metric distinguishes between feature sets that might 
have identical means but a different spread of values around the mean. In addition, Bhatta- 
charyya distance is theoretically linked to the ChernofT error bound applied to Bayesian 
classification error [Faugeras and Pratt, 1980]. If texture features are normally distributed these 
error bounds can be applied to Bhattacharyya measurements. As is true of the other metrics, 
Bhattacharyya distance is most successfully used on normalized variables. 


-61 - 


Regression Techniques 

As s first step in regression analysis, all variables were screened for asymmetrical distribu- 
tions, If histograms and quantile plots of the sorted data against quantiles of the standard normal 
distribution (Q-Q plots) were skewed, the variables were transformed using power functions 
(Tukey, 1977] that gave the best approximation to symmetric and, if possible, normal distribu- 
tions. This was done to better satisfy least-squares assumptions of normality and homoscedastic 
error. 


For each of the twelve models, the best subset of predictor variables was selected using a 
leaps and bounds regression method available from the S interactive data analysis package 
(Becker, 1984] in which Mallow's C f statistic served as 0\? criterion for goodness of fit. This 
method is a generalization of stepwise regression methods that examines all possible subsets of 
predictor variables rather than the effect of a single addition to or deletion from the predictor set. 
The C p statistic, closely related to the adjusted coefficient of determination, (Draper, 1981] is: 
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RSS f is the residual sum of squares from a model containing p parameters and a 2 is the residual 
mean square from the largest equation possible containing all the variables. 

The term a 2 is taken as the unbiased estimate of the error variance a 2 . Since the expected 
value of the C p variable is approximately p , the best, least biased equations are those in which 
the C p statistic is equal to the number of parameters. In each case, the equation with the smal- 
lest number of parameters and least biased fit was chosen for the regression model. 

Standard least-squares multiple regression was then run on the transformed variables using a 
set of programs from the S package. The significance of the regression coefficients was checked 
with a t -statistic and the overall regression significance with an F -statistic. The residual stan- 
dard error and adjusted coefficient of determination (/? 2 ) were calculated for comparison of 
models. Regression residuals were plotted against the fitted values, and a locally weighted 
smoothed line was drawn through the scatterplots to detect departures from the zero mean line. 
In addition, the sorted residuals were plotted against quartiles of the standard normal distribution 
to check for normality. Finally, a robust, iterated, weighted least-squares regression was run. 
Observations that received low weight were examined and the effect of deleting these samples was 
determined by repeating the least squares regression on the trimmed data set. If the least squares 
and robust regressions agreed well in coefficients and residuals, it was presumed that the least- 


squares assumptions had been sufficiently met and that the tests of significance stood. 

Without a large, longitudinal data set at a given site it is difficult to to test the predictive 
power of a model. In the case of the TM/MSS image pair, regression models were used for com- 
paring the effect of resolution on the scene model, not for cross-predictive purposes. Instead, for 
each model, a jackknife technique was used to verify internal stability and to test internal predic- 
tability. In this technique, each observation was deleted in turn, the regression repeated and the 
new coefficients used to predict the deleted value. The stability of each regression coefficient was 
measured by the coefficient of variation and any observations with highly deviant coefficients were 
examined for error. The predicted and obsetved values were plotted in sample sequence and as 
scatterplots. Prediction residuals were also plotted to check for any systematic patterns. The 
overall prediction performance was measured using a model efficiency statistic [Rango and Mar- 
tinec, 1079] which is a non-dimensional “goodness of fit” function: 
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Here ft is the observed snow-covered area, y is the mean snow-covered area, and ft' is the 
predicted snow-covered area. Similar to the coefficient of determination, this statistic is a meas- 
ure of the proportion of variance explained by the model. 

The purpose of independently deriving regression models for December and January images 
at a given site was to test the generality of each set of extracted parameters by means of a cross- 
prediction test. This technique determines how well one set of parameters predicts the snow 
cover at the same site but under an altered snow cover pattern. The data for the two dates were 
then pooled to obtain a general equation for the scene. As before, predictions and prediction resi- 
duals were plotted and model efficiencies calculated. 


10*0. Results 

Texture Characteristics 

Gray level co-occurrence matrices (GLCM) are joint probability tables for a specified rela- 
tionship between pixels. A comparison of matrices constrained by different joint relationships 
should reveal information about three fundamental texture properties: periodicity, directionality, 
and information content. Thus, a good starting point for analyzing the effect of resolution on 
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texture wm to compare various co-occurrence matrices and the statistics extracted from them at 
both TM and MSS resolutions. 

Using image pair A, co-occurrence matrices were tabulated over the entire TM and MSS 
scenes for a range of displacements (1, 3, 5, 10 pixels), a range of angles (0 \ 45 \ 90 \ 135 # ) and 
a range of quantisations (8, 10, 32, 04 grey levels). The chi-square statistic was used as the cri- 
terion to choose the best combination of specifications. Plots of chi-square values for all quantiza- 
tions and distances (Figure 10.0) show that the maximum chi-square angle (in this case vertical) 
was identical at all grey levels and displacements. Chi-square values from the two resolutions 
were consistently parallel in behavior; most information was found at a displacement of one pixel 
and, by definition, at the highest number of grey levels (64). 

When TM and MSS images were divided into samples (for the purpose of metric analysis) it 
was possible to compare the x 9 selected angles for each image by sample. Table 10.3 gives the 
correlations between selected angles for each image pair; only image D had a correlation better 
than 0.5. For all images, horizontal and vertical directions were dominant perhaps because the 
diagonal distance is actually 1.4 times longer than the distance to adjoining horizontal or vertical 
pixels. In the TM image the selected angle oscillated between the two orientations more rapidly 
than those of the comparable MSS image. This divergence in directionality is the most distinct 
textural difference between TM and MSS images. Although this may be interpreted simply as 
increased noise, it has elsewhere been reported that directionality is critical for distinguishing very 
similar cover types [Hallada, 1982). This implies that the tendency of investigators to reduce 
computation by averaging the various directional co-occurrence matrices is probably an unwise 
economy in classification studies. However, averaged matrices are useful for deriving rotation- 
invariant measures of texture; statistics computed from averaged matrices are the most useful 
features for predicting SCA since they hold for various orientations of the terrain image. 

Figure 10.7 is a graph of chi-square values for one angle plotted against distance. The 
exponential drop in chi-square beyond a distance of one corresponds to the sharp drop-off in the 
autocorrelation function observed in any natural texture (Laws, 1980]. Also note that large 
differences in information content due to greater quantization are only significant for distances 
less than five, beyond that, low quantization yields the same information. The large chi-square 
difference between MSS and TM images is due to the four-fold greater sample size used in the TM 
image to cover a comparable sample size at MSS resolution, not a reflection of a far greater infor- 
mation content. Regular artificial patterns will have additional peaks in the chi-square/distance 
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plot at displacement* corresponding to the periodicity of the pattern. No such underlying periodi- 
city was detected in these images using displacements that covered up to 00% of the linear 
dimension of the sample. 

By definition the \ 2 statistic increases in magnitude as the information content climbs with 
increasing quantization. The computational costs of carrying out the analysis at 64 grey levels 
were too high to be considered in this study since each doubling of the grey scale increased all cal- 
culations four-fold. Sixteen grey levels were chosen as a compromise between information content 
and efficiency based in part from the analysis of co-occurrence statistics that follows. 

GLOM Statistics 

Plots of each co occurrence statistic against angle for image pair A revealed that in each 
case the maximum chi-square angle was that which produced a matrix dominated by small grey- 
level transitions. This is a diagonally dominant matrix indicative of a relatively coarse texture 
(see Figures 10.8a,b). In short, the preferred textural orientation corresponded to the highest 
autocorrelation. 

The i ertia statistic was minimal at the x 2 selected angle because it is designed to give most 
weight to infrequent, large, high contrast, transitions. And because the grey level differences act 
to weight the statistic, most differentiation between angles was achieved at the highest quantiza- 
tion level. MSS and TM responses were parallel but MSS values were consistently higher at each 
grey level, the expected behavior for a coarser texture. 

Another statistic that varied inversely to the chi-square evaluation was entropy. This is the 
case because the negative log of small probability transitions is much greater than the negative 
log of high probability transitions. That entropy should be minimized where structure is greatest 
is intuitively correct. Entropy also dropped with decreasing quantization due to the increased 
probability per transition. Entropy values were greater for MSS than for TM at all grey levels, 
again indicating the relative coarseness of the MSS texture. 

Dominance of low “contrast” transitions was best detected by the homogeneity statistic. 
Because the gradient between joint pixels appears in the denominator of the homogeneity for- 
mula, the largest homogeneities for both TM and MSS w*re found at the lowest quantization level 
at the x 2 selected angle. Yet this statistic is relatively insensitive to quantization giving good 
differentiation between angles at the highest number of grey levels. As expected, the TM image 
had higher homogeneity than the coarser MSS image. 
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Highly corrected with homogeneity, the energy statistic wu alto maximised at lowest 
quantisation for the * 9 angle. Energy is simply the square of the transition probabilities; as the 
probability per bin drops with increasing quantisation, the square gets increasingly smaller. Thus, 
energy is very sensitive to an increase in quantisation such that no detectable difference exists 
between angles at 64 grey levels. At this quantisation the statistics would have to be scaled to 
remain compatible with other GLCM statistics. Energy wns higher for the TM scene because the 
higher probability of low contrast transitions dominated in the finer resolution image. 

Both TM and MSS images produced positive correlation (COR) statistics on the order of 
p — 0.8, an indication of a strong association between rows and columns of the GLCM. As the 
co-occurrence matrix becomes less diagonally dominant the COR value increases. Accordingly, 
the highest correlations coincided with the \ 3 angle and with the TM scene at all levels of quanti- 
sation. Unlike the energy statistic, COR, a standardized statistic ranging between zero and one, 
is useful at ail grey levels. 

Two other correlation measures based on matrix entropy are strongly associated with the 
COR statistic. Information correlation 2 (ICOR2) behaves exactly like COR although the value 
of the correlations is on average 0.1 below COR values. ICOR2 is maximal when the difference 
between total matrix entropy and the row or column entropies it, smallest i.e. when values arc 
more evenly spread throughout the matrix. Conversely, the information correlation 1 (ICORl) 
statistic assigns highest negative correlations to matrices in which this entropy differential is smal- 
lest making ICORl inversely related to COR and ICOR2. TM matrices generally have lower 
ICORl values than do MSS matrices since the difference between total and marginal entropies is 
smaller for the for the finer texture. 

It can be concluded from examination of the co-occurrence statistics that many are inter- 
correlated. Those statistics that emphasized low contrast transitions, energy and homogeneity, 
were positively correlated with each other (p = 0.82) and negatively correlated (p — -0.87) with 
entropy and inertia which give weight to high contrast transitions. The correlation statistics, 
COR, ICORl, ICOR2 were highly correlated with each other either directly or inversely. If this 
redundancy is removed, the seven co-occurrence statistics are reduced to three: measures of low 
contrast, high contrast and correlation. 

While some co-occurrence statistics arc equally effective at all levels of quantization some 
are more sensitive with fewer grey levels, others with more. A good middle ground of 10 grey lev- 
els coincides with the same choice made for the sake of computational efficiency. 
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All GLCM statistics were normalized by sample size making direct comparUor of TM and 
MSS textural features valid. The statistics in all cases are sensitive to differences in resolution 
between the scenes. MSS features consistently show a relatively coarser texture in which high 
contrast transitions form a larger proportion of concurrences than in the TM image. 

Filter Statistics 

A survey of TM and MSS Alter statistics showed that averages and sums (per unit area) are 
very close in value for the two resolutions but that the standard deviation of the TM images is 
consistently higher indicating that Altering brings forth more edge detail in the higher frequency 
TM data. In the next step of analysis, various metrics were used to quantify the resolution 
dependent differences detected by both the local statistics of the Altered images and by the co- 
occurrence statistics. 

Metric Analysis 

For each of the four sets of registered TM and MSS imagery, the distance between features 
was measured three ways. Results from using Euclidean distance, Gower Similarity, and Bhatta- 
charyya distance are given in Table 10.4. Distances between SUM feature sets are reported only 
in the case of Euclidean distance as this was the only analysis for which all variable* were normal- 
ized to a (0,1) range. The SUM variable alone is strictly dependent or. . ample size; left un~scaled, 
SUM values distorted feature space making valid Gower similarity or Bhattacharyya distance 
measures impossible. By contrast, the GLCM statistics, based on probabilities rather than fre- 
quencies and calculated from histogram equalized images arc standard statistics; likewise AVG 
and SD values are inherently normalized. A clear picture of texture separability was obtained 
despite the omission of SUM variables from Bhattacharyya distance and Gower similarity 
analysis. 

Based on mean Gower Similarity, TM and MSS textures are 91% similar when GLCM 
features arc used, 96% similar with AVG features but only 71% similar with SD features. The 
same pattern is repeated by Bhattacharyya distance measurements: SD features produce twice the 
separability of the AVG Alters, while the GLCM features fall in between. This implies, not 
surprisingly, that the ratio of within scene to between scene variance is greater for averages than 
for either standard deviations or moments of the co-occurrence matrix. 

Euclidean distances ^reed with the other two metrics by singling out the SD features as th* 
most sensitive indicators * resolution dependent textural differences. In Euclidean space, AVG 
and SUM are nearly idenu al and somewhat superior to GLCM features for texture separability. 
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Though not strictly consistent, all three metrics indicate that image pairs A and B are closer than 
images C and D. 

It can be assumed that resolution differences between TM and MSS image pairs are rela- 
tively constant since the images were taken under identical conditions. That the distances 
between pairs are not constant can be attributed to differences in registration, saturation, and sur- 
face properties between sets of images. However, if all metrics were in accord, the coefficient o» ? 
variation among the features should be fairly constant. This is the case for Euclidean distance 
and Gower similarity measures but is not for Bhattacharyya distance values, a result most likely 
due to the incorporation of the variance into the calculation of distance. 


The Bhattacharyya distance measure may also be inconsistent because it was applied to 
GLCM and filter variables that in some cases were clearly non-norm ll whereas it is defined only 
for Gaussian distributions. Likewise, the error bounds reported in Table 10.4 are not statistically 
significant but serve as rough limits on the accuracy of calculated distances. 

Metric analysis supports the conclusions gathered from inspection of individual co- 
occurrence statistics: there is a consistent, detectable difference in texture between TM and MSS 
images. This difference is on the order of 5% to 10% when characterized by GLCM statistics or 
AVG statistics but as much as 30% when the variance of the filtered images is used. Clearly, the 
difference in textural information represented by these texture features is far less than what the 
human eye perceives and what would be expected by a two-fold improvement in resolution. 
Translation of textural information into joint probabilities or into selectively filtered enhance- 
ments and subsequently into summary statistics involves a loss in information that dampens out 
distinctions between textures that are much more pronounced at the original level. It is not 
surprising that an improvement or degradation of resolution should correspond to a concomitant 
increase or decrease in image variance detectable by SD filters, but it is unexpected that the co- 
occurrence statistics are so relatively insensitive to these changes. 


Regression Analysis Results 

An initial task before undertaking regression analysis was to survey the distributions of both 
dependent and predictor variables. Summary statistics for snow-covered area (SCA) are given in 

Table 10.5. In general the standard deviations are roughly equal to — the range th’ suggesting 

4 

non-normal distributions. For a normal distribution three standard deviations on either side of 

the mean contains almost all cases making the standard deviation approximately — the range 
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(Arkin and Colton, 1070], This was confirmed by strongly skewed histograms of SCA for each 
image; low cover samples occurred far more frequently than high cover samples. This is typically 
the case for data that are counts or amounts because of a fixed zero boundary and a high or 
unlimited upper bound (Chambers, 1983]. Excessive skewness often implies a correlation of varia- 
bility with mean level which can produce non-constant variance and heteroscedastic error in a 
least-squares regression [Bartlett, 1947). Non-normality invalidates the usual significance tests 
and heteroscedaaiicity reduces the precision of the estimates. Fortunately, transformations that 
stabilize variance also tend to normalize the data; this is usually accomplished by power functions 
[Chambers, 1983]: 


*• 

e < o 

for right skew 

log y, 


6 « 0 

- Vi* 

0 > 0 

for left skew 


The best transformation was chosen by plotting sorted transformed data against quartiles of 
the standard normal distribution. The square root transformation (6 = 0.5) turned out to be best 
for the December and January data while a log transform was necessary for the MSS and TM 
image pairs. 

The predictor variables were surveyed in a similar manner for asymmetrical distributions. 
The purpose was to gain insight into the behavior of variables and to determine whether transfor- 
mation of skewed variables to symmetric improved regression models. Distributions of all the co- 
occurrence statistics show a preponderance of small transitions which is expected for a relatively 
coarse texture measured at a single pixel displacement. The inertia statistic, which compensates 
for this typical situation by giving more weight to less frequent, large transitions, is symmetric. 
Energy and homogeneity statistics are right skewed because the probability of transition term 
dominates inversely; entropy is left ikewed since the term dominates directly. By contrast, the 
response of the filter variables is scene specific; for TM and MSS images the AVG, SD and SUM 
statistics are generally symmetric while for the December and January images they are skewed 
left. 

Skewed predictor variables were transformed using power functions and then entered into 
the leaps an i bounds regression for subset selection Transformed variables do not simplify or 
improve the models. In many cases, a mix of symmetric and skewed variables have the best 
explanatory power. In final fo; n, TM/'MSS regressions are semi-log functions and the 
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December/January are “semi-square root*’ functions. 

Regression Significance 

A summary of the regression models is given in Tables 10.6a, b. The equations are grouped 
by image and type of summary statistic used for filtered samples. It is clear from these results 
that snow-covered area can be successfully regressed on a combination of texture statistics. While 
the models vary in size betw \ , ’ ree and seven parameters, the best results are obtained with 
four or five predictor variables. Ali regression coefficients are significant at the 0.005 level or 
better and the least-squares coefficients are on average within 5% of those estimated by robust 
regression. The proportion of explained variance measured by the adjusted coefficient of determi- 
nation is in all cases at least 0.95. 


Analysis of the residuals was hampered by the log transform of the dependent variable. In 
log units the residuals are homoscedastic and close to being normally distributed, once corrections 
were made for negative log values ( antilog values between zero and one). The residuals in real 
values are log-normally distributed, showing larger variance for low snow cover values. The resi- 
duals from December and January models, in square root units, were evenly spread around the 
zero mean line and close to normal on a Q-Q plot. 

Least squares linear regression assumes that regressor variables are independent, random 
variates and that errors are not autocorrelated. When regressor variables and errors are positively 
autocorrelated, the true variance of regression coefficient estimates is underestimated leading to 
overestimation of t and F significance tests and inflated R 2 values [Cliff and Ord, 1981). A basic 
property of geographic data is its spatial autocorrelation. It can be assumed, therefore, that the 
texture features used as predictors are autocorrelated to some degree. 

In order to check for residua! first-order autocorrelation Durbin-Watson statistics were cal- 
culated for all models: 


DW 


s ( - «<-! ) 2 

E u, 2 


(10.11) 


Here «q is the regression residual at location t . 

A Durbin-Watson value close to 2 indicates no autocorrelation, a value of zero implies per- 
fect positive autocorrelation and a value of 4 implies negative autocorrelation. Durbin-Watson 
statistics for all models, given in Table 10.7, are consistent: all but the December models show no 
postive first order autocorrelation (significant at the 0.01 level). Statistics for December models 




arc inconclusive; the null hypothesis of no positive autocorrelation can neither be rejected or 
accepted. Overall, it can be concluded that the sampling interval de-emphasized adjacent 
influences and failed to coincide with any large periodicities of the image function. However, in 
the general case, depending on sampling frequency and site-specific texture pattern, it may be 
necessary to include autoregressive terms for some models. 

Regression Models 

Models derived for TM and MSS images that used a combination of co-occurrence s',atistics 
and AVG filter values were very similar in composition and coefficient values. The two leading 
variables were symmetrically distributed co-occurrence statistics followed by a smoothing filter 
and a Sobel gradient operator. When AVG statistics were replaced with SUM or SL statistics, 
the models diverged in number of variables, composition and coefficient values; SD models were 
more dissimilar than AVG models. In general, the AVG/GLCM model characterizing scene tex- 
ture in the MSS image was unaltered for the TM image despite the two-fold improvement in reso- 
lution. Major differences became apparent with the use of SUM variables and were quite pro- 
nounced with SD variables. 

As can be seen from December and January imagery (Figures 10.5), the snow cover receded 
substantially in one month from 33% to 19% snow cover. Nevertheless, one pair of models 
resulted which were quite similar for both images. When GLCM and AVG features were com- 
bined, only three parameters were required and of taese two were shared in common by the 
separate models but produced very different coefficients. By contrast, the use of SUM filters 
meant a large increase in the number of parameters to seven without a corresponding improve- 
ment in R 2 . Of these seven parameters only three were in common and the coefficients were 
quite dissimilar. On the other hand, SD models for each date had five parameters, three of which 
were identical (ICORl, / 17 0 , / i/e) and two of which were highly correlated between dates. Ail 
five coefficients were quite similar but no formal test for the equality of regression coefficients 
through time could be made because the model specifications were not identical. The degree of 
similarity between models was inferred from the results of forecasting the snow cover at one date 
using model coefficients derived from the other. 

In an effort to get a ^ neral equation applicable to both December and January scenes, the 
data were pooled and re-submitted for regression analysis. The resulting models, though 
significant and predictive (see Figure 10.9) tended to be over-parameterized ranging from six pred- 
ictors using AVG filters to eight using SUM filters. Again no conclusions could be drawn about 
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the similarities of the regression relationships between pooled and individual data sets because the 
model specifications varied so widely. 

Regression Variable Selection 

Co-occurrence statistics that are designed to reflect second-order relationships did not out- 
perform the first-order filter statistics. No model was composed exclusively of co-occurrence 
statistics while the smoothing filter appeared in almost all models. Perhaps the simpler statistics 
were more suitable for use in a functional relai onship because snow-covered area is a high con- 
trast, low variance target easily captured by first-order statistics. The complicated heterogeneous 
cover of, say, an urban scene may require second-order statistics for purposes of discrimination 
and classification. 

The most significant and predictive equations were those that had four or five parameters. 
Examination of variables selected by the leaps and bounds method suggests that there are six 
categories of texture features that contain most of the texture information in a scene: low-contrast 
GLCM moments, high-contrast GLCM moments, GLCM correlations, the low pass filter (/t/ 0 ), 
vertical edge detectors and horizontal edge detectors. It is possible that a generic equation com- 
posed of variables drawn from each category but with scene specific coefficients could be generally 
applied to snow-covered watersheds. 

Prediction Results 

The significance of the regressions and the proportion of variance explained by the equations 
were all uniformly high. These models are admittedly scene-specific and calibrated with a limited 
data set. The only means of testing model predictability was to perform jackknife regressions 
summarizing internal predictability with model efficiency scores (see Table 10.8). These results 
are reported in both transformed and actual units. Clearly, the type of transformation applied to 
the dependent vari ble had a strong effect on the outcome of prediction. In log units, the 
efficiency of TM and MSS regressions was on average about 15% r than R 2 values. When 
predictions were converted to actual values (number of snow-c* ~ed pixels) the efficiencies 

dropped to zero due to excessive overshooting for the top 8% oi e values. If the domain of 
prediction is limited to low and mid-range values the efficiencies return to those measured in log 
units. The coefficients proved to be quite stable, varying on average 2.5% but in no case more 
than 6.5%. Log transformation linearized the model making possible a highly significant regres- 
sion, but because small log residuals for high values converted exponentially into much larger 
actual residuals, the domain for accurate prediction of real values was limited to 25% SCA. 
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The over-p r ediction problem inherent in log transformed data was not encountered when the 
data were square-root transformed. Jackknife efficiencies in square root units were quite high: 5% 
lower than ff 2 values and when measured in real values only 10% lower than /? 2 . Figures 
10. 10a, b and 10.1 la, b show the results of jackknife regression for the December scene in both 
actual and real values using AVG filters. 

Regressions for December and January that involved SD filteri differed from each other by 
two variables. Two edge filters in the December equation were replaced by edge/spot filters in 
the January model. The coefficients were on the whole quite similar. So it is not unexpected that 
when each model was applied to the other scene that model efficiencies, in transformed terms, 
differed by only 5%. More importantly, both models turned out to be surprisingly efficient at 
prediction: 91% and 95% of the variance was explained by the models (see Table 10.8). Figures 
10. 12a, b and 10.13a, b are plots of cross-predicted versus actual values for the two dates. 

By contrast, the AVG models for the two dates were poor predictors even though the indivi- 
dual regressions had high R 2 . The December AVG model predicted January SCA with 44% 
efficiency and the reverse was only 40% efficient. SUM models for both images were so large 
(p=7) that no attempt was made to test predictability. In general, predictive analysis re- 
emphasized what was found from jackknife regressions: standard deviations of filtered samples 
were by far the best predictors. 

The remarkably good predictability between scenes that differed in snow coverage by 15% 
implies that scene dependent parameters are robust enough to encompass recessional pattern 
changes. It may be possible to derive general parameters for moderately sized basins for use in 
the melt season when patterns of snow recession are duplicated year to year. 

10.7. Discussion and Conclusion 

Texture Characteristics 

When an image is considered to represent a random field, the cooccurrence matrix becomes 
an estimate of the joint probability density function for pixels separated by given row and column 
shifts. The autocorrelation at this spatial lag is determined by the matrix transition probabilities. 
In the images of snow cover, the joint pixel correlations at a single pixel lag were high (p=0.85) 
located along the slope of the central peak of the autocorrelation function. A comparison of co- 
occurrence statistics from different resolutions is really a matter of comparing the rate of change 
in the initial slope of the autocorrelation function; steep slopes correspond to fine textures, gentle 
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slopes to coarser textures. 

Detecting texture periodicities at far outlying joint positions is likely only for regular, 
artificial patterns. For natural scenes it is most appropriate to carry out GLCM analysis over 
relatively small windows at the minimum displacement. Judging from the results of this study, 
quantization level is most critical when using small joint displacements because it acts to enhance 
or dampen the differences in slope of the autocorrelation function. Beyond five pixels large quant- 
izations are only marginally more informative. It should also be noted that if textures are to be 
analyzed over immediately adjacent neighborhoods then directional differences between textures 
become more acute and rotationally averaged measures will overlook an important feature for dis- 
tinguishing closely related textures. 

The GLCM features are moments of the joint probability distribution that describe the 
spread of values away from the central diagonal, a feature equivalent to the autocorrelation peak 
reduced to two dimensions. As such, they are descriptors of contrast and correlation; low contrast 
transitions are close to the diagonal, high contrast transitions are more distant and the correla- 
tions reflect the degree of difference in transition values. When MSS and TM co-occurrence statis- 
tics from the same site were compared, the MSS image was coarser both in terms of contrast and 
correlation. 

Metric Analysis 

The above result only served to verify the self-evident. It was necessary to use metric 
analysis to find the magnitude and variability of these texture differences. The dissimilarity 
between the two resolutions measured between the four image pairs using the standard GLCM 
statistics or AVG/SUM statistics was only five to ten percent, but this increased to 30% when the 
standard deviation of the edge filtered samples was used instead. That the variance (second 
moment) of a local first order variable should be more informative than the mean (first moment) 
is not surprising. What is unexpected is the much greater separability afforded by a first order 
moment relative to the complex moments of the joint probability distribution. 

Laws [1980] discovered that SD measures were best for discrimination but not necessarily for 
segmentation. They acted as measures of local contrast tending to locate edges rather than 
regions. The consistently higher SD’s of the TM images no doubt indicate the greater edge detail 
created by improved resolution. As measures of local contrast SD features encompass the variabil- 
ity of an entire window whereas GLCM measures of contrast are compressed histograms of grey- 
tone transitions over a distance of a single pixel. The co-occurrence matrix provides information 
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on the adjacent (or the joint) while the filter SD provides information on the local. Gagalowicz 
[1081] and Julesz [1075] have found that texture discrimination is achieved locally, not globally so 
that properties averaged over a local region may be more informative than those ascertained from 
joint relationships. In other words, more information on texture structure may be lost by generat- 
ing the GLCM and its moments than in filtering and calculating variance. 

TM v» MSS Regressions 

Metric analysis made it clear that the relationship between resolution and texture separabil- 
ity is not one to one. In fact, it is likely that studies employing standard GLCM statistics or 
SUM/AVG filters could well substitute MSS data for TM without considerable loss in information. 
In keeping with distance metric results, the AVG and to a lesser extent SUM regression models at 
the two resolutions agreed while SD models differed considerably. When regressions were worked 
out for two dates at TM resolution, both AVG and SD models were similar in composition and, 
for the latter, in coefficient values. Any conclusions that might be drawn about resolution depen- 
dent textural differences from comparing the two sets of regression models are confounded by a 
lack of control on the SCA variable. Snov' cover varied 30% between the resolutions due to limi- 
tations of thresholding and spectral/spatial disparity between resolutions, while snow receded only 
14% between sequential images of the same resolution. 

Judging from jackknife regression results, the relationship between SCA and texture meas- 
ures was of comparable strength at the two resolutions. TM model efficiencies were only 5% 
better than MSS efficiencies. But without actual ground data to both calibrate and test the 
models it can only be said that they were both internally consistent. 

Regression and Prediction Results 

In general regression of snow area onto parameters describing its distribution was successful. 
All twelve models were highly significant and explained a large proportion of the variance. For 
December models, these results should be viewed cautiously in light of positive spatial autocorre- 
lation effects leading to underestimation of coefficient variance. Nevertheless, some of these 
models were highly predictive. At 90 - 95% efficiency the semi-square root models using SD vari- 
ables were especially good predictors. Operationally, they could be used to obtain snow cover 
estimates when snow patterns are very discontinuous making manual or digital snow classification 
difficult. 


The domain of accurate prediction is limited by the transformation chosen for the depen- 
dent variable. Transformations that stabilize variance and sy mmetricize distributions should be 
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Bdcctcd with the scale of the predicted variable in mind. For these models a square root 
transform waa preferable to taking the log. In real rather than transformed values, log models 
began to seriously overshoot at 25% snow cover while square root models were stable up to 70 or 
80% SCA. Using linear modeling and a moderate sample size this inherent problem will always 
limit the domain of prediction to some degree because the uneven distribution of sample points 
over the range of prediction results in greater uncertainty for under-represented points. 

The aim of this study was not to successfully predict SCA per ae but to establish the metho- 
dologies for eventually predicting snow depth, a more direct indicator of snowmelt runofT. To 
that end, it was learned that parameters of the texture regression models are scene and pattern 
specific but are fortunately robust enough to accommodate considerable latitude in the actual 
snow pattern formations, making them useful through time as the snowline recedes. 

It may be possible, after further studies of models derived at a variety of sites, to construct 
a standard model based on all six or some set of the six caiegories of texture features used here. 
Without question the smoothing filter ( /*7 0 ) would be a necessary parameter. Among the twelve 
models studied it was nearly ubiquitous. Correlations between /i/ 0 statistics and SCA ranged 
between 0.72 and 0.95 (see Table 10.9) making it the single most predictive variable. Obviously, 
local summary statistics of the smoothed image are highly correlated with coverage because snow 
is such a singularly bright, high-contrast target within the scene. In most models this low fre- 
quency information was augmented by high pass filter variables and some mix of GLCM statistics 
that contributed information on contrast and correlation in the immediate neighborhood of each 
pixel. In gene al, it was found that GLCM and texture energy features are well suited for detect- 
ing bright, high-contrast snow cover patterns because they operate in the first case by measuring 
contrast differences or in the latter by detecting edge structure. 

The high R 2 values of the regression models were not good indicators of predictive power. 
Using a cross-prediction test only models based on a majority of SD variables were useful under 
different conditions. Thus the results of metric and regression analysis were in accord: separabil- 
ity coincided with predictability. The standard deviation was superior for both purposes because 
it reflected the level of high frequency information consistently both within and between scenes. 
Depending on the site, AVG and SUM statistics could be highly intercorrelated or quite unrelated 
while SD values remained at a consistent level of intercorrelation at different sites. It should be 
noted that filtered images retained full quantization while GLCM images were necessarily 
reduced. The power of SD statistics may in part be due to the greater information inherent in 
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higher quantizations. It appears that breaking down texture structure into joint relationships for 
a limited number of grey tones involves a loss in information greater than that involved in taking 
local statistics of filtered images. 

Conclusion 

The difference in resolution between MSS and TM images which is so dramatic to the eye 
was critical for classification accuracy. Snow classifications from TM images were far more 
detailed and complete than MSS classifications yet principal component analysis revealed that 
spectral information did not play a large part in this disparity. The data plane used for thres- 
holding snow cover in the TM images was essentially a reduced-variance visible band somewhat 
augmented by shortwave IR data for cloud discrimination purposes. TM near infrared data were 
redundant to visible data and could be eliminated from the data plane. The MSS data plane was 
also in effect a reduced-variance visible image supplemented by near IR data which, as in TM 
component analysis, supplied mostly redundant information with respect to the first axis of max- 
imum variation. In short, for binary classifications, snow is easily detected in the visible range 
but may require shortwave IR data to aid in distinguishing cloud cover from snow. 

Measures of joint grey-tone relationships and local statistics of high frequency enhanced 
images do not duplicate perceptual sensitivity to textural detail apparent with greater resolution. 
It is not surprising that these methods do not parallel those used in complex, non-linear, hierarchi- 
cal visual processing. If the texture measures failed to achieve the sensitivity of human percep- 
tual capabilities they successfully served as descriptors of mesoscale spatial distribution in func- 
tional relationships between snow-covered area and areal distribution. Despite inherent data 
problems of skewed distributions and autocorrelated samples, linear relationships making use of 
the standard deviation of convolved images were 91- 96% efficient in predicting snow-covered 
area for a given site under two snow cover patterns. 

Exploratory and particular rather than generalized, the regression analyses presented here 
were undertaken as a feasibility study yet they did yield some general observations on the nature 
of natural textures. Texture models are a subset of image models which Ahuja and Schachter 
[1981] have grouped into the stochastic and pixel-based or the deterministic and region-based. 
Global two-dimensional stochastic models specified by a particular random field can be described 
variously by the autocorrelation function, variograms, means, gradients or spatial dependencies 
[Ahuja and Schachter, 1981). Specified in this manner, a global model can be viewed as a combi- 
nation of the ideal data modified by a point spread function plus additive noise. For example, 
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Faugeras and Pratt [1080] successfully decomposed global textures into measures of the autocorre- 
lation function and moments of the histogram of the decorrelated white noise field. This descrip- 
tion is not unique in all cases; global models breakdown when the mean and autocorrelation func- 
tion are no longer stationary. It becomes necessary to include a set of means and a spatial func- 
tion to modify the symmetric autocorrelation function. In other words, with inhomogeneities glo- 
bal models become local models. 

If global properties do not hold or if an image model cannot be formulated, s texture must 
be described empirically at the pixel or local level. Determining joint and conditional probabili- 
ties is one such approach but as Anuja and Schachter warn, a joint probability density function 
may be an overspecification lacking in abstraction. Since neither over-generalized ;!obal models, 
or highly detailed joint relationships are entirely successful descriptions of natural texture it 
becomes uecessary to describe image statistics as local spatial averages. 

On the basis of perceptual experiments, Gagalowicz [1981] found that where local second 
order statistics differ from global ones, the eye is able to detect a textural difference; visual 
discrimination is thus a local process. In addition, Julesz [1975] has speculated after years of per 
ceptual testing that visual discrimination may require only local first-order statistics of simple, 
pooled feature extractors. The results of analyzing snow cover texture support the notion that 
texture is a local property. Local standard deviations were more effective for separating textures 
and more reliable for linear prediction. Measures of local variance appear to be quite informative 
yet general enough to avoid over-sensitivity to noise. Moreover, local statistics taken from a 
series of edge enhanced images do a better job at capturing edge structure than most stochastic 
models which have been criticized for failing to account for real-world spatial structure [Modes- 
tino et al., 1981). In conclusion, the results of this study confirm those recently reported else- 
where [Pietikainen et al. t 1983) that Laws texture features are powerful and efficient descriptors of 
natural textures. 

10.8. Appendix 

Co-occurrence Features 


p (i ,/ ) = matrix entry of the normalized co-occurrence matrix. 
Ng = number of grey levels. 
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Laplaci&n Filter 

10.0. Symbol* 

AVG Average of the absolute value of pixels within a specified window in an image con- 
volved with a Laws filter. 

COR Correlation co-occurrence statistic 

ENG Energy co-occurrence statistic 

ENT Entropy co-occurrence statistic 

/ 1/ 0 Laws convolution filter 0 : smoothing filter 

/ •/ 1 Laws convolution filter 1 : vertical Sobel filter 

f il 2 Laws convolution filter 2 

/ •/ 3 Laws convolution filter 3 : horizontal Sobel filter 

/ 1 / 4 Laws convolution filter 4 

/ il b Laws convolution filter 5 

/ il 0 Laws convolution filter 6 

/ 1/ 7 Laws convolution filter 7 

f it B Laws convolution filter 8 : Laplacian filter 

GLCM Grey-level co-occurrence matrix 
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HOM Homogeneity co-occurrence statistic 

ICORl Information correlation 1 co-occunence statistic 

ICOR2 Information correlation 2 co-occurrence statistic 

INR Inertia co-occurrence statistic 

MSS Multispectral Scanner 

SCA Snow-covered area 

SD Standard deviation of pixels within a specified window of an image convolved with a 

Laws filter. 

SUM Sum of the absolute value of pixels within a specified window of an image convolved 
with a Laws filter. 

TM Thematic Mapper 




10.10. Tables 


Table 10 1 Snow Classification Training Site Statistics, DN values 


Band 

Mean 

Std Dev 

cv 

Min 

Max 

Range 

TM2 

117 

12 

10% 

42 

140 

98 

TM2/TM3 

58 

2 

4% 

53 

70 

17 

TM2/TM5 

91 

9 

10% 

34 

112 

78 

PCI (TM2,5) 

128 

11 

9% 

100 

220 

120 

PCI (TM2,3,4,5,7) 

130 

15 

12% 

101 

226 

125 

Loadings of Principal Components 


Band 

Principal Component 


PCI 

Pet Var 

PC2 Pet ’ ' 

TM2 

-052 

98 2 % 

-0 46 ! 3% 

TM3 

-0 66 


-032 

TM4 

-0 53 


082 

TM5 

-0 10 


0 14 

TM7 

-005 


0 01 

MSS4 

-047 

98 3 % 

-0 55 1.3% 

MSS5 

-0 54 


-039 

MSS6 

-0 60 


046 

MSS7 

-0 35 


057 

TM2 

-0 88 

99 1 % 

-0 15 0.9% 

TM5 

-0 15 


088 

MSS4 

-080 

97 5% 

-0.59 2 5% 

MSS7 

-0 59 


0 80 
















Bhattacharyya Distance 



Table 10 5 Snow-Covered Area 
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Tabic 10 6a Regression Models 


Image 

Filter 

Coefficients 

t value 

n 2 

F value 

DF 

Resid 

SE 

MSS 

AVG 

INR 

0 1 13 

5 58 

946 

280 

4,60 

285 



ICORl 

-5030 

-8 36 







f«o 

0011 

1087 







fil* 

-0019 

-9 39 





M r 

SD 

HOM 

1 760 

238 

955 

316 

4,60 

269 



ENT 

0 668 

354 







/«'/ 0 

0034 

9 47 







f<l t 

-0019 

-3 88 





MSS 

SUM 

INR 

5 79e-2 

375 

963 

310 

5,59 

244 



ICOR2 

3400 

6 80 







fil o 

8 53e-5 

1368 







fill 

-4 33e-5 

-3 79 







fit 3 

-4 68e-5 

-3 90 





TM 

AVG 

INR 

0 299 

7 05 

962 

410 

4,60 

315 



ICORl 

-6 86 

-891 







filo 

0017 

11 65 







fils 

-0030 

-9.98 





TM 

SD 

filo 

0064 

15.73 

982 

834 

4,60 

223 



fil 2 

-0 030 

-2.47 







filt 

-0027 

-2 16 







fils 

0051 

675 





TM 

SUM 

INR 

0.258 

601 

.965 

•142 

4,60 

304 



ICORl 

-3909 

-5 42 







filo 

4 48e-5 

13 17 







fil 4 

-2 34e-5 

-8.08 






Table 10 6b Regression Models Continued 


Image Filter Coefficients t value /? a F value DF Resid 

SE 


DEC 

AVG 

ENG 

•7 25 



COR 

-0 484 

DEC 

SD 

ICORl 

-26 88 



/•Vo 

0 102 

DEC 

SUM 

COR 

45 03 



HOM 

-32 44 

JAN 

AVG 

ENG 

44 26 



HOM 

-21 86 

JAN 

SD 

ICORl 

-1926 



f« o 

0 117 

JAN 

SUM 

COR 

16 59 



HOM 

-60 06 

DI V 7 

St 

COR 

-12 19 

JAN 


ENT 

247 


Table 10 7 Durbin-Watson Statistics 


-291 

981 

1057 

3,61 

2 165 

-2 75 
-491 

991 

1352 

5,59 

1 787 

349 
8 75 

986 

582 

7,57 

2 295 

-12 10 
703 

.950 

410 

3,61 

3 141 

-791 

-545 

988 

1036 

5,59 

1 560 

4 32 

5 80 

972 

969 

7,57 

2 484 

-985 
-3 89 

987 

1446 

7,121 

1 836 


2 47 


Image 

AVG 

SD 

SUM 


DW 

P 

DW 

P 

DW 

P 

TM 

1 66 

0 17 

2 19 

-0 09 

1 76 

0 12 

MSS 

1 61 

020 

1 66 

0 17 

1 86 

0 07 

DEC 

L39* 

030 

1,61* 

0 20 

1.60* 

020 

JAN 

1 64 

0 18 

1 68 

0 16 

2 00 

-0 002 


* cannot reject or accept null hypothesis of no positive autocorrelation. All other values show no positive 


autocorrelation at the \% level. 





. 88 • 


ORIGINAL pal :: l 
OF POOR QUALITV 


10.11. Flgur« 

Figure 10.1 Snow lexture study sites: TM image C (upper left), TM Image A (upper right), MSS 
image C (lower left), MSS Image A (lower right). 







Figure 10.3 First principal component of TM bands 2, 3, 4, 5 and 7 (upper left), first principal 
component of TM bands 2 and 5 (upper right), snow-covered area (lower left) and TM band 
2 (lower right). 



. 01 . 


igure 10.4 MSS band 4, image A (upper left), TM band 2, image A (upper right), MSS snow- 
covered area (lower left) and TM snow-covered area (lower right) 
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Figures 10.8a, b TM and MSS co-occurrence matrices, 16 grey-levels. 
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